//	Programmer	-	Sujal Bista
//	File Name	-	dataTypes.h
//	Description	:	
//					basic data types
//
//**********************************************************************

#ifndef _DATA_TYPES_H_
#define _DATA_TYPES_H_

//compile only once
#pragma once
#include <math.h>
#include <iostream>


template <typename T> class data2;
template <typename T> class data3;
template <typename T> class data4;

template <typename T> class rayData;
template <typename T> class box;
template <typename T> class quaternion;

template <typename T>
class data2
{
	public:

		union
		{
			struct
			{
				T x;
				T y;
			};

			struct
			{
				T r;
				T g;
			};

			struct
			{
				T s;
				T t;
			};

			T xy[2];
			T rg[2];
			T st[2];
		};

		//default constructor
		data2(void){}

		//constructor
		data2(const T& xx,const T& yy)
		{
			x=xx;
			y=yy;
		}

		//copy constructor
		data2(const data2& data)
		{
			x=data.x;
			y=data.y;
		}

		//copy constructor
		data2(const data3<T>& data)
		{
			x=data.x;
			y=data.y;
		}

		//copy constructor
		data2(const data4<T>& data)
		{
			x=data.x;
			y=data.y;
		}
		
		//indexing
		inline T& operator[](const unsigned int& i)
		{
			return  *((&x)+i);
		}

		//clears the data
		void clear(void)
		{
			x=y=0;
		}

		//assignment operator
		const data2& operator=(const data2& data)
		{
			x=data.x;
			y=data.y;
			return *this;
		}

		//assignment operator
		const data2& operator=(const data3<T>& data)
		{
			x=data.x;
			y=data.y;
			return *this;
		}

		//assignment operator
		const data2& operator=(const data4<T>& data)
		{
			x=data.x;
			y=data.y;
			return *this;
		}

		// data equality
		const bool operator==(const data2& vec) const
		{
			return ((x == vec.x) && (y == vec.y));
		}

		// data inequality
		const bool operator!=(const data2 &vec) const
		{
			return !(*this == vec);
		}

		// data2 add
		const data2 operator+(const data2 &vec) const
		{
			return data2(x + vec.x, y + vec.y);
		}

		// data2 add (opposite of negation)
		const data2 operator+() const
		{    
			return data2(*this);
		}

		// data2 subtraction
		const data2 operator-(const data2& vec) const
		{    
			return data2(x - vec.x, y - vec.y);
		}
     
		// data2 negation
		const data2 operator-() const
		{    
			return data2(-x, -y);
		}

		// data2 decrement
		const data2 &operator-=(const data2& vec)
		{
			x -= vec.x;
			y -= vec.y;

			return *this;
		}

		// scalar self-multiply
		const data2 &operator*=(const T &s)
		{
			x *= s;
			y *= s;
	          
			return *this;
		}

		// scalar self-divecide
		const data2 &operator/=(const T &s)
		{
			const T recip = (T)(1.0/s); // for speed, one divecision

			x *= recip;
			y *= recip;
			return *this;
		}



		// post multiply by scalar
		const data2 operator*(const T &s) const
		{
			return data2(x*s, y*s);
		}
	  
		const data2 operator*(const data2& vec) const
		{
			return data2(x*vec.x, y*vec.y);
		}


		// divide by scalar
		const data2 operator/(T s) const
		{
			s = 1/s;
			return data2(s*x, s*y);
		}

		//multiplication by T operator
		data2 operator*(const T& data)
		{
			return data2(x*data,y*data);
		}


		//self addition operator
		void operator+=(const data2& data)
		{
			x+=data.x;
			y+=data.y;
		}

		//converts vector to unit length
		void unit(void)
		{
			T multiplyer=(T)(1.0/sqrt(x*x+y*y));
			x*=multiplyer;
			y*=multiplyer;
		}

		//rotate camera along X
		void rotate(const T& theta)
		{
			T newX = (T)(cos(theta)*x - sin(theta)*y);
			T newY = (T)(sin(theta)*x + cos(theta)*y);
			x=newX;	
			y=newY;
		}


		// dot product
		friend const T dot(const data2<T> &vec1, const data2<T> &vec2)
		{
			return vec1.x*vec2.x + vec1.y*vec2.y;
		}
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const data2& vec)
		{
			out<<"<"<<vec.x<<", "<<vec.y<<">";
			return out;
		}
		//reads from the input stream
		friend std::istream& operator>>(std::istream &in, data2& vec)
		{
			in>>vec.x>>vec.y;
			return in;
		}
};
typedef data2<unsigned char> uchar2;
typedef data2<char> char2;
typedef data2<int> int2;
typedef data2<unsigned int> uint2;
typedef data2<float> float2;
typedef data2<double> double2;

template <typename T>
class data3
{
	public:

		union
		{
			struct
			{
				T x;
				T y;
				T z;
			};

			struct
			{
				T r;
				T g;
				T b;
			};

			T xyz[3];
			T rgb[3];
		};

		//default constructor
		data3(void){}

		//constructor
		data3(const T& xx,const T& yy,const T& zz)
		{
			x=xx;
			y=yy;
			z=zz;
		}

		//copy constructor
		data3(const data2<T>& data)
		{
			x=data.x;
			y=data.y;
			z=0.0;
		}

		//copy constructor
		data3(const data3& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
		}

		//copy constructor
		data3(const data4<T>& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
		}

		//indexing
		inline T& operator[](const unsigned int& i)
		{
			return  *((&x)+i);
		}

		//clear the data
		void clear(void)
		{
			x=y=z=0;
		}

		//assignment operator
		const data3& operator=(const data2<T>& data)
		{
			x=data.x;
			y=data.y;
			z=(T)(0.0);
			return *this;
		}

		//assignment operator
		const data3& operator=(const data3& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			return *this;
		}

		//assignment operator
		const data3& operator=(const data4<T>& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			return *this;
		}

		// data equality
		const bool operator==(const data3& vec) const
		{
			return ((x == vec.x) && (y == vec.y) && (z == vec.z));
		}

		// data inequality
		const bool operator!=(const data3 &vec) const
		{
			return !(*this == vec);
		}

		// data3 add
		const data3 operator+(const data3 &vec) const
		{
			return data3(x + vec.x, y + vec.y, z + vec.z);
		}

		// data3 add (opposite of negation)
		const data3 operator+() const
		{    
			return data3(*this);
		}

		// data3 subtraction
		const data3 operator-(const data3& vec) const
		{    
			return data3(x - vec.x, y - vec.y, z - vec.z);
		}
     
		// data3 negation
		const data3 operator-() const
		{    
			return data3(-x, -y, -z);
		}

		// data3 decrement
		const data3 &operator-=(const data3& vec)
		{
			x -= vec.x;
			y -= vec.y;
			z -= vec.z;

			return *this;
		}

		// scalar self-multiply
		const data3 &operator*=(const T &s)
		{
			x *= s;
			y *= s;
			z *= s;
	          
			return *this;
		}

		// scalar self-divecide
		const data3 &operator/=(const T &s)
		{
			const T recip = (T)(1.0/s); // for speed, one divecision

			x *= recip;
			y *= recip;
			z *= recip;
			return *this;
		}



		// post multiply by scalar
		const data3 operator*(const T &s) const
		{
			return data3(x*s, y*s, z*s);
		}
	  
		const data3 operator*(const data3& vec) const
		{
			return data3(x*vec.x, y*vec.y, z*vec.z);
		}


		// divide by scalar
		const data3 operator/(T s) const
		{
			s = 1/s;
			return data3(s*x, s*y, s*z);
		}

		//multiplication by T operator
		data3 operator*(const T& data)
		{
			return data3(x*data,y*data,z*data);
		}


		//self addition operator
		void operator+=(const data3& data)
		{
			x+=data.x;
			y+=data.y;
			z+=data.z;
		}

		//converts vector to unit length
		void unit(void)
		{
			T multiplyer=(T)(1.0/sqrt(x*x+y*y+z*z));
			x*=multiplyer;
			y*=multiplyer;
			z*=multiplyer;
		}

		//rotate camera along X
		void rotateX(const T& theta)
		{
			T newY = (T)(cos(theta)*y - sin(theta)*z);
			T newZ = (T)(sin(theta)*y + cos(theta)*z);
			y=newY;
			z=newZ;	
		}

		//rotate camera along Y
		void rotateY(const T& theta)
		{
			T newX = (T)(cos(theta)*x + sin(theta)*z);
			T newZ = (T)(-sin(theta)*x + cos(theta)*z);

			x=newX;
			z=newZ;
		}

		//rotate camera along Z
		void rotateZ(const T& theta)
		{
			T newX = (T)(cos(theta)*x - sin(theta)*y);
			T newY = (T)(sin(theta)*x + cos(theta)*y);
			x=newX;
			y=newY;
		}

		//borrowed from Prof. Jacobs code
		void rotate(T theta,data3 axis) 
		{
			//normalize it
			axis.unit();		

			T S[3][3] = {{   0, -axis.z,  axis.y},
							{ axis.z,    0, -axis.x},
							{-axis.y,  axis.x,    0}
							};

			T uuT[3][3] = {{axis.x*axis.x, axis.x*axis.y, axis.x*axis.z},
							{axis.y*axis.x, axis.y*axis.y, axis.y*axis.z},
							{axis.z*axis.x, axis.z*axis.y, axis.z*axis.z}};

			T newX = (T)((uuT[0][0] + cos(theta)*(1-uuT[0][0]) + sin(theta)*S[0][0])*x +
						(uuT[0][1] + cos(theta)*(0-uuT[0][1]) + sin(theta)*S[0][1])*y +
						(uuT[0][2] + cos(theta)*(0-uuT[0][2]) + sin(theta)*S[0][2])*z);

			T newY = (T)((uuT[1][0] + cos(theta)*(0-uuT[1][0]) + sin(theta)*S[1][0])*x +
						(uuT[1][1] + cos(theta)*(1-uuT[1][1]) + sin(theta)*S[1][1])*y +
						(uuT[1][2] + cos(theta)*(0-uuT[1][2]) + sin(theta)*S[1][2])*z);

			T newZ = (T)((uuT[2][0] + cos(theta)*(0-uuT[2][0]) + sin(theta)*S[2][0])*x +
						(uuT[2][1] + cos(theta)*(0-uuT[2][1]) + sin(theta)*S[2][1])*y +
						(uuT[2][2] + cos(theta)*(1-uuT[2][2]) + sin(theta)*S[2][2])*z);

			x=newX;
			y=newY;
			z=newZ;
		}

		//checks if intersects with the box
		bool doesIntersect(const box<T> &bData) const 
		{
			return bData.minX<=x && x<bData.maxX && bData.minY<=y && y<bData.maxY && bData.minZ<=z && z<bData.maxZ;
		}


		///////////////////////////////////////////////////////////////////////////////////////////////
		//computes length of vector
		friend T length(const data3& data)
		{
			return (T)sqrt(data.x*data.x+data.y*data.y+data.z*data.z);
		}


		//normalize the data
		friend data3 normalize(const data3& data)
		{
			T multiplyer=(T)(1.0/length(data));
			return data3(data.x*multiplyer,data.y*multiplyer,data.z*multiplyer);
		}


		// cross product
		friend const data3 cross(const data3 &vec1, const data3 &vec2)
		{
			return data3(vec1.y*vec2.z - vec1.z*vec2.y, vec1.z*vec2.x - vec1.x*vec2.z, vec1.x*vec2.y - vec1.y*vec2.x);
		}

		// projects vector vec1 to vec2
		friend const data3 project(const data3 &vec1, const data3 &vec2)
		{
			return (vec2*(dot(vec1,vec2)/dot(vec2,vec2)));
		}

		// dot product
		friend const T dot(const data3 &vec1, const data3 &vec2)
		{
			return vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
		}

		//reflect ray
		friend data3 reflect(const data3  &I, const data3 &N)
		{
		  return I - N * dot(N, I) * (T)(2.0);
		}

		//computes linear interpolation
		friend data3 lerp(const data3 &vec0,const data3 &vec1, const T&w)
		{
			return vec0 + (vec1-vec0)*w;
		}		

		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const data3& vec)
		{
			out<<"<"<<vec.x<<", "<<vec.y<<", "<<vec.z<<">";
			return out;
		}
		//reads from the input stream
		friend std::istream& operator>>(std::istream &in, data3& vec)
		{
			in>>vec.x>>vec.y>>vec.z;
			return in;
		}
};
typedef data3<unsigned char> uchar3;
typedef data3<char> char3;
typedef data3<int> int3;
typedef data3<unsigned int> uint3;
typedef data3<float> float3;
typedef data3<double> double3;

template <typename T>
class data4
{
	public:

		union
		{
			struct
			{
				T x;
				T y;
				T z;
				T w;
			};

			struct
			{
				T r;
				T g;
				T b;
				T a;

			};

			T xyzw[4];
			T rgba[4];
		};

		//default constructor
		data4(void){}

		//constructor
		data4(const T& xx,const T& yy,const T& zz,const T& ww)
		{
			x=xx;
			y=yy;
			z=zz;
			w=ww;
		}

		//copy constructor
		data4(const data2<T>& data)
		{
			x=data.x;
			y=data.y;
			z=(T)(0.0);
			w=(T)(0.0);
		}

		//copy constructor
		data4(const data3<T>& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			w=(T)(0.0);
		}

		//copy constructor
		data4(const data4& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			w=data.w;
		}

		//indexing
		inline T& operator[](const unsigned int& i)
		{
			return  *((&x)+i);
		}

		//clear the data
		void clear(void)
		{
			x=y=z=w=0;
		}

		//assignment operator
		const data4& operator=(const data2<T>& data)
		{
			x=data.x;
			y=data.y;
			z=(T)(0.0);
			w=(T)(0.0);
			return *this;
		}

		//assignment operator
		const data4& operator=(const data3<T>& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			w=(T)(0.0);
			return *this;
		}

		//assignment operator
		const data4& operator=(const data4& data)
		{
			x=data.x;
			y=data.y;
			z=data.z;
			w=data.w;
			return *this;
		}

		

		// data equality
		const bool operator==(const data4& vec) const
		{
			return ((x == vec.x) && (y == vec.y) && (z == vec.z) && (w == vec.w));
		}

		// data inequality
		const bool operator!=(const data4 &vec) const
		{
			return !(*this == vec);
		}

		// data4 add
		const data4 operator+(const data4 &vec) const
		{
			return data4(x + vec.x, y + vec.y, z + vec.z, w + vec.w);
		}

		// data4 add (opposite of negation)
		const data4 operator+() const
		{    
			return data4(*this);
		}

		// data4 subtraction
		const data4 operator-(const data4& vec) const
		{    
			return data4(x - vec.x, y - vec.y, z - vec.z , w - vec.w);
		}
     
		// data4 negation
		const data4 operator-() const
		{    
			return data4(-x, -y, -z, -w);
		}

		// data4 decrement
		const data4 &operator-=(const data4& vec)
		{
			x -= vec.x;
			y -= vec.y;
			z -= vec.z;
			w -= vec.w;

			return *this;
		}

		// scalar self-multiply
		const data4 &operator*=(const T &s)
		{
			x *= s;
			y *= s;
			z *= s;
			w *= s;
	          
			return *this;
		}

		// scalar self-divecide
		const data4 &operator/=(const T &s)
		{
			const T recip = (T)(1.0/s); // for speed, one divecision

			x *= recip;
			y *= recip;
			z *= recip;
			w *= recip;
			return *this;
		}



		// post multiply by scalar
		const data4 operator*(const T &s) const
		{
			return data4(x*s, y*s, z*s, w*s);
		}
	  
		const data4 operator*(const data4& vec) const
		{
			return data4(x*vec.x, y*vec.y, z*vec.z ,w*vec.w);
		}


		// divide by scalar
		const data4 operator/(T s) const
		{
			s = 1/s;
			return data4(s*x, s*y, s*z, s*w);
		}

		////multiplication by T operator
		//data4 operator*(const T& data)
		//{
		//	return data4(x*data,y*data,z*data,w*data);
		//}


		//self addition operator
		void operator+=(const data4& data)
		{
			x+=data.x;
			y+=data.y;
			z+=data.z;
			w+=data.w;
		}

		//converts vector to unit length
		void unit(void)
		{
			T multiplyer=(T)(1.0/sqrt(x*x+y*y+z*z+w*w));
			x*=multiplyer;
			y*=multiplyer;
			z*=multiplyer;
			w*=multiplyer;
		}

		//rotate camera along X
		void rotateX(const T& theta)
		{
			T newY = (T)(cos(theta)*y - sin(theta)*z);
			T newZ = (T)(sin(theta)*y + cos(theta)*z);
			y=newY;
			z=newZ;	
		}

		//rotate camera along Y
		void rotateY(const T& theta)
		{
			T newX = (T)(cos(theta)*x + sin(theta)*z);
			T newZ = (T)(-sin(theta)*x + cos(theta)*z);

			x=newX;
			z=newZ;
		}

		//rotate camera along Z
		void rotateZ(const T& theta)
		{
			T newX = (T)(cos(theta)*x - sin(theta)*y);
			T newY = (T)(sin(theta)*x + cos(theta)*y);
			x=newX;
			y=newY;
		}

		//borrowed from Prof. Jacobs code
		void rotate(T theta,data4 axis) 
		{
			//normalize it
			axis.unit();		

			T S[3][3] = {{   0, -axis.z,  axis.y},
							{ axis.z,    0, -axis.x},
							{-axis.y,  axis.x,    0}
							};

			T uuT[3][3] = {{axis.x*axis.x, axis.x*axis.y, axis.x*axis.z},
							{axis.y*axis.x, axis.y*axis.y, axis.y*axis.z},
							{axis.z*axis.x, axis.z*axis.y, axis.z*axis.z}};

			T newX = (T)((uuT[0][0] + cos(theta)*(1-uuT[0][0]) + sin(theta)*S[0][0])*x +
						(uuT[0][1] + cos(theta)*(0-uuT[0][1]) + sin(theta)*S[0][1])*y +
						(uuT[0][2] + cos(theta)*(0-uuT[0][2]) + sin(theta)*S[0][2])*z);

			T newY = (T)((uuT[1][0] + cos(theta)*(0-uuT[1][0]) + sin(theta)*S[1][0])*x +
						(uuT[1][1] + cos(theta)*(1-uuT[1][1]) + sin(theta)*S[1][1])*y +
						(uuT[1][2] + cos(theta)*(0-uuT[1][2]) + sin(theta)*S[1][2])*z);

			T newZ = (T)((uuT[2][0] + cos(theta)*(0-uuT[2][0]) + sin(theta)*S[2][0])*x +
						(uuT[2][1] + cos(theta)*(0-uuT[2][1]) + sin(theta)*S[2][1])*y +
						(uuT[2][2] + cos(theta)*(1-uuT[2][2]) + sin(theta)*S[2][2])*z);

			x=newX;
			y=newY;
			z=newZ;
		}
		
		// dot product
		friend const T dot(const data4<T> &vec1, const data4<T> &vec2)
		{
			return vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z+vec1.w*vec2.w;
		}
		//computes linear interpolation
		friend data4 lerp(const data4 &vec0,const data4 &vec1, const T&w)
		{
			return vec0 + (vec1-vec0)*w;
		}		

		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const data4& vec)
		{
			out<<"<"<<vec.x<<", "<<vec.y<<", "<<vec.z<<", "<<vec.w<<">";
			return out;
		}
		//reads from the input stream
		friend std::istream& operator>>(std::istream &in, data4& vec)
		{
			in>>vec.x>>vec.y>>vec.z>>vec.w;
			return in;
		}
};
typedef data4<unsigned char> uchar4;
typedef data4<char> char4;
typedef data4<int> int4;
typedef data4<unsigned int> uint4;
typedef data4<float> float4;
typedef data4<double> double4;

template <typename T>
class data4x4
{		
	public:
		union
		{
			struct
			{
				T        _00, _01, _02, _03;
				T        _10, _11, _12, _13;
				T        _20, _21, _22, _23;
				T        _30, _31, _32, _33;
			};
			T elements[16];
		};
				
		//default constructor
		data4x4(void){}
		
		//default constructor
		data4x4(
			const T &e0,const T &e1,const T &e2,const T &e3,
			const T &e4,const T &e5,const T &e6,const T &e7,
			const T &e8,const T &e9,const T &e10,const T &e11,
			const T &e12,const T &e13,const T &e14,const T &e15
			)
		{
			elements[0]=e0;	elements[1]=e1;	elements[2]=e2;	elements[3]=e3;
			elements[4]=e4;	elements[5]=e5;	elements[6]=e6;	elements[7]=e7;
			elements[8]=e8;	elements[9]=e9;	elements[10]=e10;	elements[11]=e11;
			elements[12]=e12;	elements[13]=e13;	elements[14]=e14;	elements[15]=e15;			
		}

		//copy const
		data4x4(const data4x4<T> &mat)
		{	
			elements[0]=mat.elements[0];
			elements[1]=mat.elements[1];
			elements[2]=mat.elements[2];
			elements[3]=mat.elements[3];

			elements[4]=mat.elements[4];
			elements[5]=mat.elements[5];
			elements[6]=mat.elements[6];
			elements[7]=mat.elements[7];

			elements[8]=mat.elements[8];
			elements[9]=mat.elements[9];
			elements[10]=mat.elements[10];
			elements[11]=mat.elements[11];

			elements[12]=mat.elements[12];
			elements[13]=mat.elements[13];
			elements[14]=mat.elements[14];
			elements[15]=mat.elements[15];
		}

		//copy const
		data4x4(const data3<T> &lookAt,const data3<T> &up,const data3<T> &right,const data3<T> &position)
		{	
			getVector(0)=right;
			getVector(1)=up;
			getVector(2)=lookAt;
			elements[3]=(T)-position.x;	
			elements[7]=(T)-position.y;
			elements[11]=(T)-position.z;
			elements[12]=(T)0.0;	elements[13]=(T)0.0;	elements[14]=(T)0.0;	elements[15]=(T)1.0;
		}


		//sets matrix to identity;
		void identity(void)
		{
			elements[0]=(T)1.0;	elements[1]=(T)0.0;	elements[2]=(T)0.0;	elements[3]=(T)0.0;
			elements[4]=(T)0.0;	elements[5]=(T)1.0;	elements[6]=(T)0.0;	elements[7]=(T)0.0;
			elements[8]=(T)0.0;	elements[9]=(T)0.0;	elements[10]=(T)1.0;	elements[11]=(T)0.0;
			elements[12]=(T)0.0;	elements[13]=(T)0.0;	elements[14]=(T)0.0;	elements[15]=(T)1.0;
		}

		//default constructor
		void clear(void)
		{
			memeset(elements,0,sizeof(T)*16);
		}

		//returns right
		inline void setRight(const data3<T> &right)
		{
			getVector(0)=right;
		}

		//returns up
		inline void setUp(const data3<T> &up)
		{
			getVector(1)=up;
		}

		//returns lookat
		inline void setLookAt(const data3<T> &lookAt)
		{
			getVector(2)=lookAt;
		}		

		//returns position
		inline void setPosition(const data3<T> &position)
		{
			elements[3]=(T)-position.x;	
			elements[7]=(T)-position.y;
			elements[11]=(T)-position.z;
		}

		//returns right
		inline const data3<T> &getRight(void)
		{
			return getVector(0);
		}

		//returns up
		inline const data3<T> &getUp(void)
		{
			return getVector(1);
		}

			//returns lookat
		inline const data3<T> &getLookAt(void)
		{
			return getVector(2);
		}		

		//returns position
		inline const data3<T> getPosition(void)
		{
			return data3<T>(elements[3],elements[7],elements[11]);
		}

		//indexing
		inline const T& operator[](const unsigned int& i) const
		{
			return elements[i];
		}

		//transposes a matrix
		inline data4x4<T> transpose(void) const
		{
			data4x4<T> dataTranspose;
			dataTranspose.elements[0]=elements[0];
			dataTranspose.elements[1]=elements[4];
			dataTranspose.elements[2]=elements[8];
			dataTranspose.elements[3]=elements[12];
			dataTranspose.elements[4]=elements[1];
			dataTranspose.elements[5]=elements[5];
			dataTranspose.elements[6]=elements[9];
			dataTranspose.elements[7]=elements[13];
			dataTranspose.elements[8]=elements[2];
			dataTranspose.elements[9]=elements[6];
			dataTranspose.elements[10]=elements[10];
			dataTranspose.elements[11]=elements[14];
			dataTranspose.elements[12]=elements[3];
			dataTranspose.elements[13]=elements[7];
			dataTranspose.elements[14]=elements[11];
			dataTranspose.elements[15]=elements[15];
			return dataTranspose;
		}

		//inverse of the Rotation translatio matrix
		inline data4x4<T> inverseRT(void) const
		{
			data4x4<T> invTranslation(
				1, 0, 0, -elements[3],
				0, 1, 0, -elements[7],
				0, 0, 1, -elements[11],
				0, 0, 0, 1
				);

			data4x4<T> invRotation(
				elements[0],	elements[4],	elements[8],	0,
				elements[1],	elements[5],	elements[9],	0,
				elements[2],	elements[6],	elements[10],0,
				0,	0,	0,	elements[15]
			);
			return invTranslation*invRotation;
		}

		//compute inverse of general matrix
		inline data4x4<T> inverse(void) const
		{
			data4x4<T> inverseMatrix;
			//compute the determinant
			T det=determinant();

			inverseMatrix.elements[0]= (-elements[13]*elements[10]*elements[7] +elements[9]*elements[14]*elements[7] +elements[13]*elements[6]*elements[11]-elements[5]*elements[14]*elements[11] -elements[9]*elements[6]*elements[15] +elements[5]*elements[10]*elements[15])/det;
			inverseMatrix.elements[4]= ( elements[12]*elements[10]*elements[7] -elements[8]*elements[14]*elements[7] -elements[12]*elements[6]*elements[11]+elements[4]*elements[14]*elements[11] +elements[8]*elements[6]*elements[15] -elements[4]*elements[10]*elements[15])/det;
			inverseMatrix.elements[8]= (-elements[12]*elements[9]* elements[7] +elements[8]*elements[13]*elements[7] +elements[12]*elements[5]*elements[11]-elements[4]*elements[13]*elements[11] -elements[8]*elements[5]*elements[15] +elements[4]*elements[9]* elements[15])/det;
			inverseMatrix.elements[12]=( elements[12]*elements[9]* elements[6] -elements[8]*elements[13]*elements[6] -elements[12]*elements[5]*elements[10]+elements[4]*elements[13]*elements[10] +elements[8]*elements[5]*elements[14] -elements[4]*elements[9]* elements[14])/det;
			
			inverseMatrix.elements[1]= ( elements[13]*elements[10]*elements[3] -elements[9]*elements[14]*elements[3] -elements[13]*elements[2]*elements[11]+elements[1]*elements[14]*elements[11] +elements[9]*elements[2]*elements[15] -elements[1]*elements[10]*elements[15])/det;
			inverseMatrix.elements[5]= (-elements[12]*elements[10]*elements[3] +elements[8]*elements[14]*elements[3] +elements[12]*elements[2]*elements[11]-elements[0]*elements[14]*elements[11] -elements[8]*elements[2]*elements[15] +elements[0]*elements[10]*elements[15])/det;
			inverseMatrix.elements[9]= ( elements[12]*elements[9]* elements[3] -elements[8]*elements[13]*elements[3] -elements[12]*elements[1]*elements[11]+elements[0]*elements[13]*elements[11] +elements[8]*elements[1]*elements[15] -elements[0]*elements[9]* elements[15])/det;
			inverseMatrix.elements[13]=(-elements[12]*elements[9]* elements[2] +elements[8]*elements[13]*elements[2] +elements[12]*elements[1]*elements[10]-elements[0]*elements[13]*elements[10] -elements[8]*elements[1]*elements[14] +elements[0]*elements[9]* elements[14])/det;
			
			inverseMatrix.elements[2]= (-elements[13]*elements[6]* elements[3] +elements[5]*elements[14]*elements[3] +elements[13]*elements[2]*elements[7]-elements[1]*elements[14]*elements[7] -elements[5]*elements[2]*elements[15] +elements[1]*elements[6]* elements[15])/det;
			inverseMatrix.elements[6]= ( elements[12]*elements[6]* elements[3] -elements[4]*elements[14]*elements[3] -elements[12]*elements[2]*elements[7]+elements[0]*elements[14]*elements[7] +elements[4]*elements[2]*elements[15] -elements[0]*elements[6]* elements[15])/det;
			inverseMatrix.elements[10]=(-elements[12]*elements[5]* elements[3] +elements[4]*elements[13]*elements[3] +elements[12]*elements[1]*elements[7]-elements[0]*elements[13]*elements[7] -elements[4]*elements[1]*elements[15] +elements[0]*elements[5]* elements[15])/det;
			inverseMatrix.elements[14]=( elements[12]*elements[5]* elements[2] -elements[4]*elements[13]*elements[2] -elements[12]*elements[1]*elements[6]+elements[0]*elements[13]*elements[6] +elements[4]*elements[1]*elements[14] -elements[0]*elements[5]* elements[14])/det;
			
			inverseMatrix.elements[3]= ( elements[9]* elements[6]* elements[3] -elements[5]*elements[10]*elements[3] -elements[9]* elements[2]*elements[7]+elements[1]*elements[10]*elements[7] +elements[5]*elements[2]*elements[11] -elements[1]*elements[6]* elements[11])/det;
			inverseMatrix.elements[7]= (-elements[8]* elements[6]* elements[3] +elements[4]*elements[10]*elements[3] +elements[8]* elements[2]*elements[7]-elements[0]*elements[10]*elements[7] -elements[4]*elements[2]*elements[11] +elements[0]*elements[6]* elements[11])/det;
			inverseMatrix.elements[11]=( elements[8]* elements[5]* elements[3] -elements[4]*elements[9]* elements[3] -elements[8]* elements[1]*elements[7]+elements[0]*elements[9]* elements[7] +elements[4]*elements[1]*elements[11] -elements[0]*elements[5]* elements[11])/det;
			inverseMatrix.elements[15]=(-elements[8]* elements[5]* elements[2] +elements[4]*elements[9]* elements[2] +elements[8]* elements[1]*elements[6]-elements[0]*elements[9]* elements[6] -elements[4]*elements[1]*elements[10] +elements[0]*elements[5]* elements[10])/det;

			return inverseMatrix;

		}

				
		//returns vectors
		inline data3<T> & getVector(const int& i)
		{
			return *((data3<T> *)&elements[(i*4)]);
		}

		//multiply
		inline data3<T> operator*(const data3<T> &vec) const
		{
			return data3<T>(
				elements[0]*vec.x+elements[1]*vec.y+elements[2]*vec.z,
				elements[4]*vec.x+elements[5]*vec.y+elements[6]*vec.z,
				elements[8]*vec.x+elements[9]*vec.y+elements[10]*vec.z				
				);
		}

		//multiply
		inline data4<T> operator*(const data4<T> &vec) const
		{
			return data4<T>(
				elements[0]*vec.x+elements[1]*vec.y+elements[2]*vec.z+elements[3]*vec.w,
				elements[4]*vec.x+elements[5]*vec.y+elements[6]*vec.z+elements[7]*vec.w,
				elements[8]*vec.x+elements[9]*vec.y+elements[10]*vec.z+elements[11]*vec.w,
				elements[12]*vec.x+elements[13]*vec.y+elements[14]*vec.z+elements[15]*vec.w
				);
		}

		//multiply
		inline data4x4<T> operator*(const data4x4<T> &mat) const
		{
			data4x4<T> data(
			
				elements[0]*mat.elements[0]+elements[1]*mat.elements[4]+elements[2]*mat.elements[8]+elements[3]*mat.elements[12],
				elements[0]*mat.elements[1]+elements[1]*mat.elements[5]+elements[2]*mat.elements[9]+elements[3]*mat.elements[13],
				elements[0]*mat.elements[2]+elements[1]*mat.elements[6]+elements[2]*mat.elements[10]+elements[3]*mat.elements[14],
				elements[0]*mat.elements[3]+elements[1]*mat.elements[7]+elements[2]*mat.elements[11]+elements[3]*mat.elements[15],

				elements[4]*mat.elements[0]+elements[5]*mat.elements[4]+elements[6]*mat.elements[8]+elements[7]*mat.elements[12],
				elements[4]*mat.elements[1]+elements[5]*mat.elements[5]+elements[6]*mat.elements[9]+elements[7]*mat.elements[13],
				elements[4]*mat.elements[2]+elements[5]*mat.elements[6]+elements[6]*mat.elements[10]+elements[7]*mat.elements[14],
				elements[4]*mat.elements[3]+elements[5]*mat.elements[7]+elements[6]*mat.elements[11]+elements[7]*mat.elements[15],

				elements[8]*mat.elements[0]+elements[9]*mat.elements[4]+elements[10]*mat.elements[8]+elements[11]*mat.elements[12],
				elements[8]*mat.elements[1]+elements[9]*mat.elements[5]+elements[10]*mat.elements[9]+elements[11]*mat.elements[13],
				elements[8]*mat.elements[2]+elements[9]*mat.elements[6]+elements[10]*mat.elements[10]+elements[11]*mat.elements[14],
				elements[8]*mat.elements[3]+elements[9]*mat.elements[7]+elements[10]*mat.elements[11]+elements[11]*mat.elements[15],

				elements[12]*mat.elements[0]+elements[13]*mat.elements[4]+elements[14]*mat.elements[8]+elements[15]*mat.elements[12],
				elements[12]*mat.elements[1]+elements[13]*mat.elements[5]+elements[14]*mat.elements[9]+elements[15]*mat.elements[13],
				elements[12]*mat.elements[2]+elements[13]*mat.elements[6]+elements[14]*mat.elements[10]+elements[15]*mat.elements[14],
				elements[12]*mat.elements[3]+elements[13]*mat.elements[7]+elements[14]*mat.elements[11]+elements[15]*mat.elements[15]
			);
			
			return data;
		}

		//multiply
		inline void operator*=(const data4x4<T> &mat)
		{	

			float4 temp=float4(elements[0],elements[1],elements[2],elements[3]);
			elements[0]=temp.x*mat.elements[0]+temp.y*mat.elements[4]+temp.z*mat.elements[8]+temp.w*mat.elements[12];
			elements[1]=temp.x*mat.elements[1]+temp.y*mat.elements[5]+temp.z*mat.elements[9]+temp.w*mat.elements[13];
			elements[2]=temp.x*mat.elements[2]+temp.y*mat.elements[6]+temp.z*mat.elements[10]+temp.w*mat.elements[14];
			elements[3]=temp.x*mat.elements[3]+temp.y*mat.elements[7]+temp.z*mat.elements[11]+temp.w*mat.elements[15];

			temp=float4(elements[4],elements[5],elements[6],elements[7]);
			elements[4]=temp.x*mat.elements[0]+temp.y*mat.elements[4]+temp.z*mat.elements[8]+temp.w*mat.elements[12];
			elements[5]=temp.x*mat.elements[1]+temp.y*mat.elements[5]+temp.z*mat.elements[9]+temp.w*mat.elements[13];
			elements[6]=temp.x*mat.elements[2]+temp.y*mat.elements[6]+temp.z*mat.elements[10]+temp.w*mat.elements[14];
			elements[7]=temp.x*mat.elements[3]+temp.y*mat.elements[7]+temp.z*mat.elements[11]+temp.w*mat.elements[15];

			temp=float4(elements[8],elements[9],elements[10],elements[11]);
			elements[8]=temp.x*mat.elements[0]+temp.y*mat.elements[4]+temp.z*mat.elements[8]+temp.w*mat.elements[12];
			elements[9]=temp.x*mat.elements[1]+temp.y*mat.elements[5]+temp.z*mat.elements[9]+temp.w*mat.elements[13];
			elements[10]=temp.x*mat.elements[2]+temp.y*mat.elements[6]+temp.z*mat.elements[10]+temp.w*mat.elements[14];
			elements[11]=temp.x*mat.elements[3]+temp.y*mat.elements[7]+temp.z*mat.elements[11]+temp.w*mat.elements[15];

			temp=float4(elements[12],elements[13],elements[14],elements[15]);
			elements[12]=temp.x*mat.elements[0]+temp.y*mat.elements[4]+temp.z*mat.elements[8]+temp.w*mat.elements[12];
			elements[13]=temp.x*mat.elements[1]+temp.y*mat.elements[5]+temp.z*mat.elements[9]+temp.w*mat.elements[13];
			elements[14]=temp.x*mat.elements[2]+temp.y*mat.elements[6]+temp.z*mat.elements[10]+temp.w*mat.elements[14];
			elements[15]=temp.x*mat.elements[3]+temp.y*mat.elements[7]+temp.z*mat.elements[11]+temp.w*mat.elements[15];			
		}

		//translate the matrix
		inline const data4x4<T> &translate(const data3<T> &data)
		{
			//make the translate matrix
			data4x4<T> dataMat(			
				1, 0, 0, data.x,
				0, 1, 0, data.y,
				0, 0, 1, data.z,
				0, 0, 0, 1
				);

			//multiply the matrix
			(*this)*=(dataMat);

			//return
			return (*this);
		}

		//rotate camera along X
		inline const data4x4<T> &rotateX(const T& theta)
		{
			//make the translate matrix
			data4x4<T> dataMat(			
				1, 0, 0, 0,
				0, (T)cos(theta),(T)-sin(theta), 0,
				0, (T)sin(theta),(T)cos(theta), 0,
				0, 0, 0, 1
				);

			//multiply the matrix
			(*this)*=(dataMat);

			//return
			return (*this);
		}

		//rotate camera along Y
		inline const data4x4<T> &rotateY(const T& theta)
		{
			//make the translate matrix
			data4x4<T> dataMat(			
				(T)cos(theta), 0,(T)sin(theta), 0,
				0, 1, 0, 0,
				(T)-sin(theta), 0,(T)cos(theta), 0,
				0, 0, 0, 1
				);

			//multiply the matrix
			(*this)*=(dataMat);

			//return
			return (*this);
		}

		//rotate camera along Z
		inline const data4x4<T> &rotateZ(const T& theta)
		{
			//make the translate matrix
			data4x4<T> dataMat(			
				(T)cos(theta),-(T)sin(theta), 0, 0,
				(T)sin(theta), (T)cos(theta), 0, 0,
				0, 0, 1, 0,
				0, 0, 0, 1
				);

			//multiply the matrix
			(*this)*=(dataMat);

			//return
			return (*this);
		}

		//scale the matrix
		inline const data4x4<T> &scale(const data3<T> &data)
		{
			//make the translate matrix
			data4x4<T> dataMat(			
				data.x, 0, 0, 0,
				0, data.y, 0, 0,
				0, 0, data.z, 0,
				0, 0, 0, 1
				);

			//multiply the matrix
			(*this)*=(dataMat);

			//return
			return (*this);
		}

		inline T determinant(void) const
		{
			T retVal=
					elements[0] * elements[5] * elements[10] * elements[15] + elements[0] * elements[6] * elements[11] * elements[13] + elements[0] * elements[7] * elements[9] * elements[14]
					+ elements[1] * elements[4] * elements[11] * elements[14] + elements[1] * elements[6] * elements[8] * elements[15] + elements[1] * elements[7] * elements[10] * elements[12]
					+ elements[2] * elements[4] * elements[9] * elements[15] + elements[2] * elements[5] * elements[11] * elements[12] + elements[2] * elements[7] * elements[8] * elements[13]
					+ elements[3] * elements[4] * elements[10] * elements[13] + elements[3] * elements[5] * elements[8] * elements[14] + elements[3] * elements[6] * elements[9] * elements[12]
					- elements[0] * elements[5] * elements[11] * elements[14] - elements[0] * elements[6] * elements[9] * elements[15] - elements[0] * elements[7] * elements[10] * elements[13]
					- elements[1] * elements[4] * elements[10] * elements[15] - elements[1] * elements[6] * elements[11] * elements[12] - elements[1] * elements[7] * elements[8] * elements[14]
					- elements[2] * elements[4] * elements[11] * elements[13] - elements[2] * elements[5] * elements[8] * elements[15] - elements[2] * elements[7] * elements[9] * elements[12]
					- elements[3] * elements[4] * elements[9] * elements[14] - elements[3] * elements[5] * elements[10] * elements[12] - elements[3] * elements[6] * elements[8] * elements[13];

			return retVal;
		}
		//gets the rotaion part of the matrix
		const data4x4<T> &operator=(const quaternion<T> &Q)
		{
			//From Quaternion to Matrix and Back
			//J.M.P. van Waveren
			//ID software
			//http://cache-www.intel.com/cd/00/00/29/37/293748_293748.pdf

			const T *q = (T*) &Q.x;
			T x2 = q[0] + q[0];
			T y2 = q[1] + q[1];
			T z2 = q[2] + q[2];
			T w2 = q[3] + q[3];
			T yy2 = q[1] * y2;
			T xy2 = q[0] * y2;
			T xz2 = q[0] * z2;
			T yz2 = q[1] * z2;
			T zz2 = q[2] * z2;
			T wz2 = q[3] * z2;
			T wy2 = q[3] * y2;
			T wx2 = q[3] * x2;
			T xx2 = q[0] * x2;
		
			elements[0]=(T)(- yy2 - zz2 + 1.0);		elements[1]=xy2 + wz2;					elements[2]=xz2 - wy2;					elements[3]=(T)0.0;
			elements[4]=xy2 - wz2;					elements[5]=(T)(- xx2 - zz2 + 1.0);		elements[6]=yz2 + wx2;					elements[7]=(T)0.0;
			elements[8]=xz2 + wy2;					elements[9]=yz2 - wx2;					elements[10]=(T)(-xx2 - yy2 + 1.0);		elements[11]=(T)0.0;
			elements[12]=(T)0.0;					elements[13]=(T)0.0;					elements[14]=(T)0.0;					elements[15]=(T)1.0;

			return *this;
		}


		//transposes a matrix
		friend data4x4<T> transpose(const data4x4<T>& data)
		{
			return data.transpose();
		}
		//extracts euler angles from rotation matrix
		//0 = x,y,z
		//1 = x,z,y
		//2 = y,x,z
		//3 = y,z,x
		//4 = z,x,y
		//5 = z,y,x
		friend data3<T> extractEulerAngles(data4x4<T> &rotationMatrix,unsigned int order)
		{
			data3<T> theta;
			if(order==0) //x,y,z
			{
				//////////////////////////////////////////////////////////////////////////
				//borrowed from http://www.geometrictools.com/Documentation/EulerAngles.pdf
				//converts the rotation matrix into xyz rotation order
				if (rotationMatrix.elements[0*4+2] < +1)
				{
					if (rotationMatrix.elements[0*4+2] > -1)
					{
						theta.y = (T)asin(rotationMatrix.elements[0*4+2]);
						theta.x = (T)atan2(-rotationMatrix.elements[1*4+2],rotationMatrix.elements[2*4+2]);
						theta.z = (T)atan2(-rotationMatrix.elements[0*4+1],rotationMatrix.elements[0*4+0]);
					}
					else // rotationMatrix.elements[2] = -1
					{
						// Not a unique solution: theta.z - theta.x = atan2(rotationMatrix.elements[4],rotationMatrix.elements[5])
						theta.y = (T)-PI/2;
						theta.x = (T)-atan2(rotationMatrix.elements[1*4+0],rotationMatrix.elements[1*4+1]);
						theta.z = (T)0;
					}
				}
				else // rotationMatrix.elements[2] = +1
				{
					// Not a unique solution: theta.z + theta.x = atan2(rotationMatrix.elements[4],rotationMatrix.elements[5])
					theta.y = (T)+PI/2;
					theta.x = (T)atan2(rotationMatrix.elements[1*4+0],rotationMatrix.elements[1*4+1]);
					theta.z = (T)0;
				}
			}
			else if(order==1) //xzy
			{
				if (rotationMatrix.elements[1] < +1)
				{
					if (rotationMatrix.elements[1] > -1)
					{
						theta.z = (T)asin(-rotationMatrix.elements[1]);
						theta.x = (T)atan2(rotationMatrix.elements[9],rotationMatrix.elements[5]);
						theta.y = (T)atan2(rotationMatrix.elements[2],rotationMatrix.elements[0]);
					}
					else // rotationMatrix.elements[1] = -1
					{
						// Not a unique solution: theta.y - theta.x = atan2(-rotationMatrix.elements[8],rotationMatrix.elements[10])
						theta.z = (T)+PI/2;
						theta.x = (T)atan2(-rotationMatrix.elements[8],rotationMatrix.elements[10]);
						theta.y = (T)0;
					}
				}
				else // rotationMatrix.elements[1] = +1
				{
					// Not a unique solution: theta.y + theta.x = atan2(-rotationMatrix.elements[8],rotationMatrix.elements[10])
					theta.z = (T)-PI/2;
					theta.x = (T)atan2(-rotationMatrix.elements[8],rotationMatrix.elements[10]);
					theta.y = (T)0;
				}
			}
			else if(order==2)	//y,x,z
			{
				if (rotationMatrix.elements[6] < +1)
				{
					if (rotationMatrix.elements[6] > -1)
					{
						theta.x = (T)asin(-rotationMatrix.elements[6]);
						theta.y = (T)atan2(rotationMatrix.elements[2],rotationMatrix.elements[10]);
						theta.z = (T)atan2(rotationMatrix.elements[4],rotationMatrix.elements[5]);
					}
					else // rotationMatrix.elements[6] = -1
					{
						// Not a unique solution: theta.z - theta.y = atan2(-rotationMatrix.elements[1],rotationMatrix.elements[0])
						theta.x = (T)+PI/2;
						theta.y = (T)-atan2(-rotationMatrix.elements[1],rotationMatrix.elements[0]);
						theta.z = (T)0;
					}
				}
				else // rotationMatrix.elements[6] = +1
				{
					// Not a unique solution: theta.z + theta.y = atan2(-rotationMatrix.elements[1],rotationMatrix.elements[0])
					theta.x = (T)-PI/2;
					theta.y = (T)atan2(-rotationMatrix.elements[1],rotationMatrix.elements[0]);
					theta.z = (T)0;
				}
			}
			else if(order==3) //y,z,x
			{
				if (rotationMatrix.elements[4] < +1)
				{
					if (rotationMatrix.elements[4] > -1)
					{
						theta.z = (T)asin(rotationMatrix.elements[4]);
						theta.y = (T)atan2(-rotationMatrix.elements[8],rotationMatrix.elements[0]);
						theta.x = (T)atan2(-rotationMatrix.elements[6],rotationMatrix.elements[5]);
					}
					else // rotationMatrix.elements[4] = -1
					{
						// Not a unique solution: theta.x - theta.y = atan2(rotationMatrix.elements[9],rotationMatrix.elements[10])
						theta.z = (T)-PI/2;
						theta.y = (T)-atan2(rotationMatrix.elements[9],rotationMatrix.elements[10]);
						theta.x = (T)0;
					}
				}
				else
				{
					// Not a unique solution: theta.x + theta.y = atan2(rotationMatrix.elements[9],rotationMatrix.elements[10])
					theta.z = (T)+PI/2;
					theta.y = (T)atan2(rotationMatrix.elements[9],rotationMatrix.elements[10]);
					theta.x = (T)0;
				}
			}
			else if(order==4) //z,x,y
			{

				if (rotationMatrix.elements[9] < +1)
				{
					if (rotationMatrix.elements[9] > -1)
					{
						theta.x = (T)asin(rotationMatrix.elements[9]);
						theta.z = (T)atan2(-rotationMatrix.elements[1],rotationMatrix.elements[5]);
						theta.y = (T)atan2(-rotationMatrix.elements[8],rotationMatrix.elements[10]);
					}
					else // rotationMatrix.elements[9] = -1
					{
						// Not a unique solution: theta.y - theta.z = atan2(rotationMatrix.elements[2],rotationMatrix.elements[0])
						theta.x = (T)-PI/2;
						theta.z = (T)-atan2(rotationMatrix.elements[2],rotationMatrix.elements[0]);
						theta.y = (T)0;
					}
				}
				else // rotationMatrix.elements[9] = +1
				{
					// Not a unique solution: theta.y + theta.z = atan2(rotationMatrix.elements[2],rotationMatrix.elements[0])
					theta.x = (T)+PI/2;
					theta.z = (T)atan2(rotationMatrix.elements[2],rotationMatrix.elements[0]);
					theta.y = (T)0;
				}
			}
			else if(order==5)//z,y,x
			{
				if (rotationMatrix.elements[8] < +1)
				{
					if (rotationMatrix.elements[8] > -1)
					{
						theta.y = (T)asin(-rotationMatrix.elements[8]);
						theta.z = (T)atan2(rotationMatrix.elements[4],rotationMatrix.elements[0]);
						theta.x = (T)atan2(rotationMatrix.elements[9],rotationMatrix.elements[10]);
					}
					else // rotationMatrix.elements[8] = -1
					{
						// Not a unique solution: theta.x - theta.z = atan2(-rotationMatrix.elements[6],rotationMatrix.elements[5])
						theta.y = (T)+PI/2;
						theta.z = (T)atan2(-rotationMatrix.elements[6],rotationMatrix.elements[5]);
						theta.x = (T)0;
					}
				}
				else // rotationMatrix.elements[8] = +1
				{
					// Not a unique solution: theta.x + theta.z = atan2(-rotationMatrix.elements[6],rotationMatrix.elements[5])
					theta.y = (T)-PI/2;
					theta.z = (T)atan2(-rotationMatrix.elements[6],rotationMatrix.elements[5]);
					theta.x = (T)0;
				}
			}
			return theta;
		}
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out,data4x4<T> &mat)
		{
			out<<mat[0]<<"\t"<<mat[1]<<"\t"<<mat[2]<<"\t"<<mat[3]<<"\t"<<"\n";
			out<<mat[4]<<"\t"<<mat[5]<<"\t"<<mat[6]<<"\t"<<mat[7]<<"\t"<<"\n";
			out<<mat[8]<<"\t"<<mat[9]<<"\t"<<mat[10]<<"\t"<<mat[11]<<"\t"<<"\n";
			out<<mat[12]<<"\t"<<mat[13]<<"\t"<<mat[14]<<"\t"<<mat[15]<<"\t"<<std::endl;
			return out;
		}
		//reads from the input stream
		friend std::istream& operator>>(std::istream &in,data4x4<T> &mat)
		{
			in>>mat.elements[0]>>mat.elements[1]>>mat.elements[2]>>mat.elements[3];
			in>>mat.elements[4]>>mat.elements[5]>>mat.elements[6]>>mat.elements[7];
			in>>mat.elements[8]>>mat.elements[9]>>mat.elements[10]>>mat.elements[11];
			in>>mat.elements[12]>>mat.elements[13]>>mat.elements[14]>>mat.elements[15];
			return in;
		}
		//translate the matrix
		friend const data4x4<T> translate(const data3<T> &data)
		{
			//make the translate matrix
			return data4x4<T>(			
				1, 0, 0, data.x,
				0, 1, 0, data.y,
				0, 0, 1, data.z,
				0, 0, 0, 1
				);
		}
		//rotate camera along X
		friend const data4x4<T> rotateX(const T& theta)
		{
			//make the translate matrix
			return data4x4<T>(			
				1, 0, 0, 0,
				0, (T)cos(theta),(T)-sin(theta), 0,
				0, (T)sin(theta),(T)cos(theta), 0,
				0, 0, 0, 1
				);
		}
		//rotate camera along Y
		friend const data4x4<T> rotateY(const T& theta)
		{
			//make the translate matrix
			return  data4x4<T>(			
				(T)cos(theta), 0,(T)sin(theta), 0,
				0, 1, 0, 0,
				(T)-sin(theta), 0,(T)cos(theta), 0,
				0, 0, 0, 1
				);
		}
		//rotate camera along Z
		friend const data4x4<T> rotateZ(const T& theta)
		{
			//make the translate matrix
			return data4x4<T>(			
				(T)cos(theta),-(T)sin(theta), 0, 0,
				(T)sin(theta), (T)cos(theta), 0, 0,
				0, 0, 1, 0,
				0, 0, 0, 1
				);
		}
		//scale the matrix
		friend const data4x4<T> scale(const data3<T> &data)
		{
			//make the translate matrix
			return data4x4<T>(			
				data.x, 0, 0, 0,
				0, data.y, 0, 0,
				0, 0, data.z, 0,
				0, 0, 0, 1
				);
		}
};

typedef data4x4<int> int4x4;
typedef data4x4<float> float4x4;
typedef data4x4<double> double4x4;

//quaternion class
template <typename T>
class quaternion
{
	union
		{
			struct
			{
				T x;
				T y;
				T z;
				T w;
			};
			struct
			{
				data3<T> axis;
				T scale;
			};
			T xyzw[4];
		};
	public:
		quaternion(){}

		quaternion(const quaternion &q)
		{	
			x=q.x;
			y=q.y;
			z=q.z;
			w=q.w;
		}

		quaternion(const T& nx,const T& ny,const T& nz,const T& nw)
		{	
			x=nx;
			y=ny;
			z=nz;
			w=nw;
		}

		quaternion(const data3<T> &v, const T &w):axis(v),scale(w)
		{
		
		}

		//clear the data
		void clear(void)
		{
			x=0;
			y=0;
			z=0;
			w=1;		
		}

		//conjugate of quaterion
		//same as inverse aslong as the length is unit
		quaternion conjugate() const
		{
			return quaternion(-x,-y,-z,w);
		}

		//converts vector to unit length
		void unit(void)
		{
			T multiplyer=(T)(1.0/sqrt(x*x+y*y+z*z+w*w));
			x*=multiplyer;
			y*=multiplyer;
			z*=multiplyer;
			w*=multiplyer;
		}

		quaternion operator*(const quaternion &q) const
		{
			//multiply the quaternion
			return quaternion(w * q.x + x * q.w + y * q.z - z * q.y,
							  w * q.y + y * q.w + z * q.x - x * q.z,
							  w * q.z + z * q.w + x * q.y - y * q.x,
							  w * q.w - x * q.x - y * q.y - z * q.z);

		}

		const quaternion &operator*=(const quaternion &q)
		{
			*this=*this*q;
			return *this;
		}

		// multiplying a quaternion q with a vector v applies the q-rotation to v
		data3<T> operator* (const data3<T> &vec) const
		{
			//create a quat
			quaternion vecQuat(normalize(vec),0);
	 
			//v=q * v *conjugate(q)
			quaternion resQuat = *this * vecQuat * conjugate();
			return data3<T>(resQuat.x, resQuat.y, resQuat.z);
		}

		// Convert from Axis Angle
		void fromAxis(const data3<T> &v, const T &angle)
		{
			axis=normalize(v)*(T)(sin(angle*0.5));
			w = (T)cos(angle*0.5);
		}

		// Convert from Euler Angles
		void fromEulerAngles(const T &pitch, const T &yaw,const T &roll)
		{

			//fixed
			//http://gpwiki.org/index.php/OpenGL:Tutorials:Using_Quaternions_to_represent_rotation

			T sinp = (T)sin(-pitch*.5);
			T siny = (T)sin(-yaw*.5);
			T sinr = (T)sin(-roll*.5);
			T cosp = (T)cos(-pitch*.5);
			T cosy = (T)cos(-yaw*.5);
			T cosr = (T)cos(-roll*.5);

			x = sinp * cosy * cosr - cosp * siny * sinr;
			y = cosp * siny * cosr + sinp * cosy * sinr;
			z = cosp * cosy * sinr - sinp * siny * cosr;
			w = cosp * cosy * cosr + sinp * siny * sinr;

			unit();
		}

		// Convert to Axis/Angles
		void getAxisAngle(data3<T> &axis, T &angle)
		{
			T scale = sqrt(x * x + y * y + z * z);
			axis.x = x / scale;
			axis.y = y / scale;
			axis.z = z / scale;
			angle = (T)(acos(w) * 2.0);
		}

		// quaternion add
		const quaternion operator+(const quaternion &vec) const
		{
			return quaternion(x + vec.x, y + vec.y, z + vec.z, w + vec.w);
		}

		// quaternion subtraction
		const quaternion operator-(const quaternion& vec) const
		{    
			return quaternion(x - vec.x, y - vec.y, z - vec.z , w - vec.w);
		}

		// add (opposite of negation)
		const quaternion operator+() const
		{    
			return quaternion(*this);
		}
 		// negation
		const quaternion operator-() const
		{    
			return quaternion(-x, -y, -z, -w);
		}

		// post multiply by scalar
		const quaternion operator*(const T &s) const
		{
			return quaternion(x*s, y*s, z*s, w*s);
		}

		T dot(const quaternion & b) const
		{
			return x*b.x+y*b.y+z*b.z+w*b.w;
		}

		//gets the rotaion part of the matrix
		const quaternion &operator=(data4x4<T> &matrix)
		{
			////normalize the data
			//matrix.getVector(0).unit();
			//matrix.getVector(1).unit();
			//matrix.getVector(2).unit();

			//From Quaternion to Matrix and Back
			//J.M.P. van Waveren
			//ID software
			//http://cache-www.intel.com/cd/00/00/29/37/293748_293748.pdf

			T s0, s1, s2;
			int k0, k1, k2, k3;
			T *q = (T*) &x;
			const T *m = (T*) &matrix.elements[0];

			if ( m[0 * 4 + 0] + m[1 * 4 + 1] + m[2 * 4 + 2] > (T)0.0 ) 
			{
				k0 = 3;
				k1 = 2;
				k2 = 1;
				k3 = 0;
				s0 = (T)1.0;
				s1 = (T)1.0;
				s2 = (T)1.0;
			} 
			else if ( m[0 * 4 + 0] > m[1 * 4 + 1] && m[0 * 4 + 0] > m[2 * 4 + 2] ) 
			{
				k0 = 0;
				k1 = 1;
				k2 = 2;
				k3 = 3;
				s0 = (T)1.0;
				s1 = (T)-1.0;
				s2 = (T)-1.0;
			} 
			else if (m[1 * 4 + 1] > m[2 * 4 + 2]) 
			{
				k0 = 1;
				k1 = 0;
				k2 = 3;
				k3 = 2;
				s0 = (T)-1.0;
				s1 = (T)1.0;
				s2 = (T)-1.0;
			} 
			else
			{
				k0 = 2;
				k1 = 3;
				k2 = 0;
				k3 = 1;
				s0 = (T)-1.0;
				s1 = (T)-1.0;
				s2 = (T)1.0;
			}
			T t = s0 * m[0 * 4 + 0] + s1 * m[1 * 4 + 1] + s2 * m[2 * 4 + 2] + (T)1.0;
			T s = (T)(0.5/sqrt(t));
			q[k0] = s * t;
			q[k1] = ( m[0 * 4 + 1] - s2 * m[1 * 4 + 0] ) * s;
			q[k2] = ( m[2 * 4 + 0] - s1 * m[0 * 4 + 2] ) * s;
			q[k3] = ( m[1 * 4 + 2] - s0 * m[2 * 4 + 1] ) * s;
			return *this;
		}

		
		// Convert to Matrix
		data4x4<T> matrix(void) const
		{
			//From Quaternion to Matrix and Back
			//J.M.P. van Waveren
			//ID software
			//http://cache-www.intel.com/cd/00/00/29/37/293748_293748.pdf

			const T *q = (T*) &x;
			T x2 = q[0] + q[0];
			T y2 = q[1] + q[1];
			T z2 = q[2] + q[2];
			//T w2 = q[3] + q[3];
			T yy2 = q[1] * y2;
			T xy2 = q[0] * y2;
			T xz2 = q[0] * z2;
			T yz2 = q[1] * z2;
			T zz2 = q[2] * z2;
			T wz2 = q[3] * z2;
			T wy2 = q[3] * y2;
			T wx2 = q[3] * x2;
			T xx2 = q[0] * x2;

			return data4x4<T>(
				(T)(- yy2 - zz2 + 1.0),		xy2 + wz2, 					xz2 - wy2, 				(T)0.0,
				xy2 - wz2,					(T)(- xx2 - zz2 + 1.0),		yz2 + wx2,				(T)0.0,
				xz2 + wy2,					yz2 - wx2,					(T)(-xx2 - yy2 + 1.0),	(T)0.0,
				(T)0.0,						(T)0.0,						(T)0.0,					(T)1.0
			);
		}

		/////////////////////////////////////////////////////////////////////////////////////////////
		//spherical linear interpolation
		// dot product
		friend const T dot(const quaternion<T> &vec1, const quaternion<T> &vec2)
		{
			return vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z+vec1.w*vec2.w;
		}
		friend quaternion<T> slerp(const quaternion<T> &v0, const quaternion<T> &v1, const T &w)
		{
			//find cos angles
			T cosAngle = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z+v0.w*v1.w;
			
			//check if the two quaternions are very close
			if (cosAngle > (T)0.99995) 
			{
				quaternion<T> result=v0+(v1-v0)*w;
				result.unit();
				return result;
			}

			//clamp the values
			cosAngle=MAX(cosAngle,-1);
			cosAngle=MIN(cosAngle,1);

			T theta = (T)acos(cosAngle); //angle between input quaternion
			T slerpTheta = (T)(theta*w);    //position of the interpolation

			quaternion<T> v2 = (v1-v0*cosAngle); //make v2  orthogonal basis on the direction of v1
			v2.unit();//normalize it

			//perform spherical interpolation
			return v0*cos(slerpTheta) + v2*sin(slerpTheta);
		}
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const quaternion& vec)
		{
			out<<"<"<<vec.x<<", "<<vec.y<<", "<<vec.z<<", "<<vec.w<<">";
			return out;
		}
		//reads from the input stream
		friend std::istream& operator>>(std::istream &in, quaternion& vec)
		{
			in>>vec.x>>vec.y>>vec.z>>vec.w;
			return in;
		}

}; 

typedef quaternion<int> intQuaternion;
typedef quaternion<float> floatQuaternion;
typedef quaternion<double> doubleQuaternion;

//line data packet
template <typename T>
class line2
{
	public:
		data2<T> points[2];

		line2(void)
		{
		}

		line2(const data2<T> &pos1, const data2<T> &pos2)
		{
			points[0]=pos1;
			points[1]=pos2;
		}

		//clears the data
		void clear(void)
		{
			points[0].clear();
			points[1].clear();
		}

		//distance from point to line
		float distance(const data2<T> &pos)
		{
			T c1,c2;

			data2<T> v=points[1]-points[0];
			data2<T> w=pos-points[0];

			//calculate the projectio of p1 onto  the line segment
			c1=dot(w,v);
			if(c1<=0)
				return length(pos-points[0]);

			c2=dot(v,v);
			if(c2<=c1)
				return length(pos-points[1]);

			//how far is the projection from c1 [0-1]
			c1=c1/c2;

			//simple paramatic eqn
			w=points[0]+v*c1;

			return length(pos-w);
		}

		//project point onto the line
		//if the pos gets projected outside the line
		//it gets clamped to the endpoint
		data2<T> projectClamped(const data2<T> &pos)
		{
			T c1,c2;

			data2<T> v=points[1]-points[0];
			data2<T> w=pos-points[0];

			//calculate the projectio of p1 onto  the line segment
			c1=dot(w,v);
			if(c1<=0)
				return points[0];

			c2=dot(v,v);
			if(c2<=c1)
				return points[1];

			//how far is the projection from c1 [0-1]
			c1=c1/c2;

			//simple paramatic eqn
			w=points[0]+v*c1;

			return (w);
		}

		////checks if the triange intersects with the ray
		//bool doesIntersect(const line2<data2<T>> &l) const
		//{
		//	T t0,t1,tI;

		//	const data2<T> &startA=points[0];
		//	data2<T> dirA=points[1]-points[0];
		//	t0=(double)length(dirA);
		//	if(t0>.00001)
		//		dirA*=(1.0/t0);


		//	const data2<T> &startB=l.points[0];
		//	data2<T> dirB=l.points[1]-l.points[0];
		//	t1=(double)length(dirB);
		//	if(t1>.00001)
		//		dirB*=(1.0/t1);

		//	

		//	if(abs(dirA.x*dirB.y-dirA.y*dirB.x)>.00001)
		//	{
		//		//(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)=(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)
		//		tI = (double)-(dirA.x*startB.y-dirA.y*startB.x-dirA.x*dirA.y*startA.y+dirA.x*dirA.y*startA.x)/(dirA.x*dirB.y-dirA.y*dirB.x);
		//	}
		//	else if(abs(dirA.y*dirB.z-dirA.z*dirB.y)>.00001)
		//	{
		//		//(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)=(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)
		//		tI = (double)-(dirA.y*startB.z-dirA.z*startB.y-dirA.y*dirA.z*startA.z+dirA.y*dirA.z*startA.y)/(dirA.y*dirB.z-dirA.z*dirB.y); 
		//	}
		//	else if(abs(dirA.x*dirB.z-dirA.z*dirB.x)>.00001)
		//	{
		//		//(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)=(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)
		//		tI = (double)-(dirA.x*startB.z-dirA.z*startB.x-dirA.x*dirA.z*startA.z+dirA.x*dirA.z*startA.x)/(dirA.x*dirB.z-dirA.z*dirB.x);
		//	}
		//	else
		//		return false;

		//	//compute intersection point
		//	data2<T> intersectionPoint=startB+dirB*tI;

		//	if(tI<=t1)
		//		return true;

		//	if(length(intersectionPoint-startA)<=t0)
		//		return true;		

		//	return true;
		//}
};
typedef line2<int> int2Line;
typedef line2<float> float2Line;
typedef line2<double> double2Line;

//line data packet
template <typename T>
class line3
{
	public:
		data3<T> points[2];

		line3(void)
		{
		}

		line3(const data3<T> &pos1, const data3<T> &pos2)
		{
			points[0]=pos1;
			points[1]=pos2;
		}

		//clears the data
		void clear(void)
		{
			points[0].clear();
			points[1].clear();
		}

		//distance from point to line
		T distance(const data3<T> &pos)
		{
			T c1,c2;

			data3<T> v=points[1]-points[0];
			data3<T> w=pos-points[0];

			//calculate the projectio of p1 onto  the line segment
			c1=dot(w,v);
			if(c1<=0)
				return length(pos-points[0]);

			c2=dot(v,v);
			if(c2<=c1)
				return length(pos-points[1]);

			//how far is the projection from c1 [0-1]
			c1=c1/c2;

			//simple paramatic eqn
			w=points[0]+v*c1;

			return (T)length(pos-w);
		}

		//returns distance from point[0] to the nearest projection
		T projectDistance(const data3<T> &pos)
		{
			T c1,c2;

			data3<T> v=points[1]-points[0];
			data3<T> w=pos-points[0];

			//calculate the projectio of p1 onto  the line segment
			c1=dot(w,v);
			c2=dot(v,v);
			c1=c1/c2;

			//return the distance
			return (c1/c2);
		}

		//project point onto the line
		//if the pos gets projected outside the line
		//it gets clamped to the endpoint
		data3<T> projectClamped(const data3<T> &pos)
		{
			T c1,c2;

			data3<T> v=points[1]-points[0];
			data3<T> w=pos-points[0];

			//calculate the projectio of p1 onto  the line segment
			c1=dot(w,v);
			if(c1<=0)
				return points[0];

			c2=dot(v,v);
			if(c2<=c1)
				return points[1];

			//how far is the projection from c1 [0-1]
			c1=c1/c2;

			//simple paramatic eqn
			w=points[0]+v*c1;

			return (w);
		}

		//checks if the triange intersects with the ray
		bool doesIntersect(const line3<T> &l) const
		{
			T t0,t1,tI;

			const data3<T> &startA=points[0];
			data3<T> dirA=points[1]-points[0];
			t0=(T)length(dirA);
			if(t0>(T).00001)
				dirA*=(T)(1.0/t0);

			const data3<T> &startB=l.points[0];
			data3<T> dirB=l.points[1]-l.points[0];
			t1=(T)length(dirB);
			if(t1>(T).00001)
				dirB*=(T)(1.0/t1);
			
			if(abs(dirA.x*dirB.y-dirA.y*dirB.x)>(T).00001)
			{
				//(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)=(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)
				tI = (T)-(dirA.x*startB.y-dirA.y*startB.x-dirA.x*dirA.y*startA.y+dirA.x*dirA.y*startA.x)/(dirA.x*dirB.y-dirA.y*dirB.x);
			}
			else if(abs(dirA.y*dirB.z-dirA.z*dirB.y)>(T).00001)
			{
				//(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)=(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)
				tI = (T)-(dirA.y*startB.z-dirA.z*startB.y-dirA.y*dirA.z*startA.z+dirA.y*dirA.z*startA.y)/(dirA.y*dirB.z-dirA.z*dirB.y); 
			}
			else if(abs(dirA.x*dirB.z-dirA.z*dirB.x)>(T).00001)
			{
				//(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)=(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)
				tI = (T)-(dirA.x*startB.z-dirA.z*startB.x-dirA.x*dirA.z*startA.z+dirA.x*dirA.z*startA.x)/(dirA.x*dirB.z-dirA.z*dirB.x);
			}
			else
				return false;
			
			if(tI>t1)
				return false;

			//compute intersection point
			data3<T> intersectionPoint=startB+dirB*tI;

			if(length(intersectionPoint-startA)>t0)
				return false;

			return true;
		}

		//checks if the triange intersects with the ray
		bool getIntersect(const line3<T> &l,data3<T> &intersectionPoint) const
		{
			T t0,t1,tI;

			const data3<T> &startA=points[0];
			data3<T> dirA=points[1]-points[0];
			t0=(T)length(dirA);
			if(t0>(T).00001)
				dirA*=(T)(1.0/t0);

			const data3<T> &startB=l.points[0];
			data3<T> dirB=l.points[1]-l.points[0];
			t1=(T)length(dirB);
			if(t1>(T).00001)
				dirB*=(T)(1.0/t1);
			
			if(abs(dirA.x*dirB.y-dirA.y*dirB.x)>(T).00001)
			{
				//(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)=(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)
				tI = (T)-(dirA.x*startB.y-dirA.y*startB.x-dirA.x*dirA.y*startA.y+dirA.x*dirA.y*startA.x)/(dirA.x*dirB.y-dirA.y*dirB.x);
			}
			else if(abs(dirA.y*dirB.z-dirA.z*dirB.y)>(T).00001)
			{
				//(startB.y+(tB*dirB.y)-(startA.y*dirA.y))/(dirA.y)=(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)
				tI = (T)-(dirA.y*startB.z-dirA.z*startB.y-dirA.y*dirA.z*startA.z+dirA.y*dirA.z*startA.y)/(dirA.y*dirB.z-dirA.z*dirB.y); 
			}
			else if(abs(dirA.x*dirB.z-dirA.z*dirB.x)>(T).00001)
			{
				//(startB.z+(tB*dirB.z)-(startA.z*dirA.z))/(dirA.z)=(startB.x+(tB*dirB.x)-(startA.x*dirA.x))/(dirA.x)
				tI = (T)-(dirA.x*startB.z-dirA.z*startB.x-dirA.x*dirA.z*startA.z+dirA.x*dirA.z*startA.x)/(dirA.x*dirB.z-dirA.z*dirB.x);
			}
			else
				return false;
			
			if(tI>t1)
				return false;

			//compute intersection point
			intersectionPoint=startB+dirB*tI;

			if(length(intersectionPoint-startA)>t0)
				return false;

			return true;
		}

		//returns length of the line
		friend T length(const line3& lData)
		{
			data3<T> data=lData.points[1]-lData.points[0];
			return (T)sqrt(data.x*data.x+data.y*data.y+data.z*data.z);
		}
};
typedef line3<int> int3Line;
typedef line3<float> float3Line;
typedef line3<double> double3Line;

//line data packet
template <typename T>
class plane
{
	public:
		data3<T> normal;
		T d;

		plane(void){}


		plane(const T &a,const T &b,const T &c,const T &d):normal(a,b,c)
		{
			this->d=d;
		}

		plane(const data3<T> &p1,const data3<T> &p2,const data3<T> &p3)
		{
			normal=cross((p2-p1),(p3-p1));
			normal.unit();
			d=-dot(p1,normal);
		}

		plane(const data3<T> &n,const T &d):normal(n)
		{
			this->d=d;
		}

		plane(const data3<T> &n,const data3<T> &p1):normal(n)
		{
			d=-dot(p1,normal);
		}

		//clears the data
		void clear(void)
		{
			d=0;
			normal.clear();
		}
				
		//distance from point to plane with abs function
		float distanceRaw(const data3<T> &pos) const
		{
			return (dot(pos,normal)+d);
		}


		//distance from point to plane
		float distance(const data3<T> &pos) const
		{
			return abs(dot(pos,normal)+d);
		}

		//check if the point is infront
		bool isFront(const data3<T> &pos) const
		{
			if(distanceRaw(pos)>=0)
				return true;
			return false;
		}

		//returns the closest point in the plane
		data3<T> closestPoint(const data3<T> &pos) const
		{
			return (pos-normal*distanceRaw(pos));
		}

		//checks if the triange intersects with the ray
		bool doesIntersect(const rayData<T> &ray) const
		{
			T t=-(dot(normal,ray.position)+d)/(dot(normal,ray.direction));
			if(t<(T)(-0.00001))
				return false;
			return true;
		}

		//get returns the intersection distance
		T getIntersectDistance(const rayData<T> &ray) const
		{
			T t=-(dot(normal,ray.position)+d)/(dot(normal,ray.direction));
			if(t<(T)(-0.00001))
				return (T)-1.0;
			return t;
		}

		// post multiply by scalar
		const plane<T> operator*(const T &s) const
		{
			return  plane<T>(normal*s, d*s);
		}
		
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const plane& vec)
		{
			out<<"<"<<vec.normal.x<<", "<<vec.normal.y<<", "<<vec.normal.z<<", "<<vec.d<<">";
			return out;
		}
};
typedef plane<int> intPlane;
typedef plane<float> floatPlane;
typedef plane<double> doublePlane;


//ray data packet
template <typename T>
class rayData
{
	public:
		data3<T> position;
		data3<T> direction;

		rayData(void)
		{
		}

		rayData(const data3<T> &position, const data3<T> &direction)
		{
			this->position=position;
			this->direction=direction;
		}

		//clears the data
		void clear(void)
		{
			position.clear();
			direction.clear();
		}
};
typedef rayData<int> intRayData;
typedef rayData<float> floatRayData;
typedef rayData<double> doubleRayData;

template <typename T>
class photonData
{
	public:
		unsigned int index;
		data3<T> position;       // position ( 3 x 32 bit floats )
		data3<T> incidentDir;	//incident direction
		data3<T> color;        //color   

		photonData(void)
		{

		}

		photonData(const data3<T> &position, const data3<T> &incidentDir, const data3<T> &color)
		{
			this->position=position;
			this->incidentDir=incidentDir;
			this->color=color;
		}

		//clears the data
		void clear(void)
		{
			position.clear();
			incidentDir.clear();
			color.clear();
		}

		void setIndex(const unsigned int &index)
		{
			this->index=index;
		}

		const data3<T> &getPosition(void)
		{
			return position;
		}

		//returns difference
		data3<T> operator-(const data3<T> &data) const
		{
			return position-data;
		}

		//returns difference
		bool operator<(const photonData<T> &data) const
		{
			return (index<data.index);
		}

};
typedef photonData<int> intPhotonData;
typedef photonData<float> floatPhotonData;
typedef photonData<double> doublePhotonData;



//forward declaration
template <typename T> class box;

template <typename T>
class triangle
{
	public:
		data3<T> vertices[3];

		triangle(void){}

		triangle(const data3<T> &vertex1,const data3<T> &vertex2,const data3<T> &vertex3)
		{
			vertices[0]=vertex1;
			vertices[1]=vertex2;
			vertices[2]=vertex3;
		}

		//clears the data
		void clear(void)
		{
			vertices[0].clear();
			vertices[1].clear()
			vertices[2].clear()
		}

		//returns the normal of the triangle
		data3<T> getNormal(void) const
		{
			//compute edge
			data3<T> edge1=vertices[1]-vertices[0];
			data3<T> edge2=vertices[2]-vertices[0];

			//compute normal
			return normalize(cross(edge1,edge2));
		}

		//smaller than operator
		bool operator<(const triangle<T> &data) const
		{
			return false;
		}

		// data equality
		const bool operator==(const triangle<T> &data) const
		{
			return (vertices[0]==data.vertices[0] && vertices[1]==data.vertices[1] && vertices[2]==data.vertices[2]);
		}
		
		int planeBoxOverlap(data3<T> &normal, data3<T> &vert, data3<T> &maxbox) const // -NJMP-
		{
			data3<T> vmin,vmax;

			if(normal.x>(T)(0.0))
			{
				vmin.x=-maxbox.x - vert.x; // -NJMP-
				vmax.x= maxbox.x - vert.x; // -NJMP-
			}
			else
			{
				vmin.x= maxbox.x - vert.x; // -NJMP-
				vmax.x=-maxbox.x - vert.x; // -NJMP-
			}

			if(normal.y>(T)(0.0))
			{
				vmin.y=-maxbox.y - vert.y; // -NJMP-
				vmax.y= maxbox.y - vert.y; // -NJMP-
			}
			else
			{
				vmin.y= maxbox.y - vert.y; // -NJMP-
				vmax.y=-maxbox.y - vert.y; // -NJMP-
			}

			if(normal.z>(T)(0.0))
			{
				vmin.z=-maxbox.z - vert.z; // -NJMP-
				vmax.z= maxbox.z - vert.z; // -NJMP-
			}
			else
			{
				vmin.z= maxbox.z - vert.z; // -NJMP-
				vmax.z=-maxbox.z - vert.z; // -NJMP-
			}

			if(dot(normal,vmin)>(T)(0.0)) 
				return false;  // -NJMP-
			if(dot(normal,vmax)>=(T)(0.0)) 
				return true; // -NJMP-

			return false;
		}
		
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const triangle& data)
		{
			out<<"{"<<data.vertices[0]<<","<<data.vertices[1]<<","<<data.vertices[2]<<"}";
			return out;
		}

		//checks if the triange intersects with the ray
		bool doesIntersect(const rayData<T> &ray) const
		{
			//compute edge
			data3<T> edge1=vertices[1]-vertices[0];
			data3<T> edge2=vertices[2]-vertices[0];

			data3<T> pvec=cross(ray.direction,edge2);
			T det=dot(edge1,pvec);

			//check if the detrminant is zero
			if(det>.00001 && det <.00001)
				return false;

			T invDet=(T)(1.0/det);
			data3<T> tvec=ray.position-vertices[0];
			T u = dot(tvec, pvec) * invDet;

			if (u < 0.0 || u > 1.0)
				return false;

			//prepare to test V parameter
			data3<T> qvec=cross(tvec, edge1);

			//calculate V parameter and test bounds
			T v = dot(ray.direction, qvec) * invDet;
			if (v < 0.0 || u + v > 1.0)
				return false;

			//calculate t, ray intersects triangle
			T t = dot(edge2, qvec) * invDet;

			if(t<(T)(-0.00001))
				return false;
			return true;
		}

		
		//returns barycentric coordinate
		data3<T> getBarycentric(const data3<T> &pos) const
		{
			// Compute the normal of the triangle
			data3<T> normal = normalize(cross(vertices[1]-vertices[0],vertices[2]-vertices[0]));

			// Compute twice area of triangle ABC
			T areaABC = (T)dot(normal,cross(vertices[1]-vertices[0],vertices[2]-vertices[0]));

			// Compute a
			T areaPBC = (T)dot(normal,cross(vertices[1]-pos,vertices[2]-pos));
			T a = (T) (areaPBC / areaABC);

			// Compute b
			T areaPCA = (T)dot(normal,cross(vertices[2]-pos,vertices[0]-pos));
			T b = (T)(areaPCA / areaABC);

			return data3<T>((T)a,(T)b,(T)(1.0-a-b));
		}

		//get returns the intersection distance
		T getIntersectDistance(const rayData<T> &ray) const
		{
			//compute edge
			data3<T> edge1=vertices[1]-vertices[0];
			data3<T> edge2=vertices[2]-vertices[0];

			//compute normal
			data3<T> normal=normalize(cross(edge1,edge2));

			data3<T> pvec=cross(ray.direction,edge2);
			T det=dot(edge1,pvec);

			//check if the detrminant is zero
			if(det>.00001 && det <.00001)
				return (T)-1;

			T invDet=(T)(1.0/det);
			data3<T> tvec=ray.position-vertices[0];
			T u = dot(tvec, pvec) * invDet;

			if (u < 0.0 || u > 1.0)
				return (T)-1;

			//prepare to test V parameter
			data3<T> qvec=cross(tvec, edge1);

			//calculate V parameter and test bounds
			T v = dot(ray.direction, qvec) * invDet;
			if (v < 0.0 || u + v > 1.0)
				return (T)-1;

			//calculate t, ray intersects triangle
			T t = dot(edge2, qvec) * invDet;

			return t;
		}



		//======================== X-tests ========================
		#define AXISTEST_X01(a, b, fa, fb)           \
		   p0 = a*v0.y - b*v0.z;                   \
		   p2 = a*v2.y - b*v2.z;                   \
		   if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
		   rad = fa * boxHalfSize.y + fb * boxHalfSize.z;   \
		   if (min>rad || max<-rad) return false;

		#define AXISTEST_X2(a, b, fa, fb)            \
		   p0 = a*v0.y - b*v0.z;                  \
		   p1 = a*v1.y - b*v1.z;                   \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.y + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		//======================== Y-tests ========================
		#define AXISTEST_Y02(a, b, fa, fb)           \
		   p0 = -a*v0.x + b*v0.z;                  \
		   p2 = -a*v2.x + b*v2.z;                        \
		   if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		#define AXISTEST_Y1(a, b, fa, fb)            \
		   p0 = -a*v0.x + b*v0.z;                  \
		   p1 = -a*v1.x + b*v1.z;                     \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		//======================== Z-tests ========================
		#define AXISTEST_Z12(a, b, fa, fb)           \
		   p1 = a*v1.x - b*v1.y;                  \
		   p2 = a*v2.x - b*v2.y;                   \
		   if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.y;   \
		   if(min>rad || max<-rad) return false;


		#define AXISTEST_Z0(a, b, fa, fb)            \
		   p0 = a*v0.x - b*v0.y;             \
		   p1 = a*v1.x - b*v1.y;                  \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.y;   \
		   if(min>rad || max<-rad) return false;

		#define FINDMINMAX(x0,x1,x2,min,max) \
		   min = max = x0;   \
		   if(x1<min) min=x1;\
		   if(x1>max) max=x1;\
		   if(x2<min) min=x2;\
		   if(x2>max) max=x2;

		bool doesIntersect(const box<T> &boxData) const
		{
			data3<T> boxHalfSize=boxData.getHalfSize();
			data3<T> boxCenter=boxData.getCenter();

			//use separating axis theorem to test overlap between triangle and box 
			//    need to test for overlap in these directions: 
			//    1) the {x,y,z}-directions (actually, since we use the AABB of the triangle 
			//       we do not even need to test these) 
			//    2) normal of the triangle 
			//    3) crossproduct(edge from tri, {x,y,z}-directin) 
			//       this gives 3x3=9 more tests 
			data3<T> v0,v1,v2;
			T min,max,p0,p1,p2,rad,fex,fey,fez;      // -NJMP- "d" local variable removed
			data3<T> normal,e0,e1,e2;

			// This is the fastest branch on Sun 
			// move everything so that the boxCenter is in (0,0,0) 
			v0=vertices[0]-boxCenter;
			v1=vertices[1]-boxCenter;
			v2=vertices[2]-boxCenter;

			// compute triangle edges 
			e0=v1-v0;      // tri edge 0 
			e1=v2-v1;      // tri edge 1 
			e2=v0-v2;      // tri edge 2 

			// Bullet 3:  
			//  test the 9 tests first (this was faster) 
			fex = abs(e0.x);
			fey = abs(e0.y);
			fez = abs(e0.z);
			AXISTEST_X01(e0.z, e0.y, fez, fey);
			AXISTEST_Y02(e0.z, e0.x, fez, fex);
			AXISTEST_Z12(e0.y, e0.x, fey, fex);

			fex = abs(e1.x);
			fey = abs(e1.y);
			fez = abs(e1.z);
			AXISTEST_X01(e1.z, e1.y, fez, fey);
			AXISTEST_Y02(e1.z, e1.x, fez, fex);
			AXISTEST_Z0(e1.y, e1.x, fey, fex);

			fex = abs(e2.x);
			fey = abs(e2.y);
			fez = abs(e2.z);
			AXISTEST_X2(e2.z, e2.y, fez, fey);
			AXISTEST_Y1(e2.z, e2.x, fez, fex);
			AXISTEST_Z12(e2.y, e2.x, fey, fex);

			// Bullet 1: 
			//  first test overlap in the {x,y,z}-directions 
			//  find min, max of the triangle each direction, and test for overlap in 
			//  that direction -- this is equivalent to testing a minimal AABB around 
			//  the triangle against the AABB 

			// test in X-direction 
			FINDMINMAX(v0.x,v1.x,v2.x,min,max);
			if(min>boxHalfSize.x || max<-boxHalfSize.x)
				return false;

			// test in Y-direction 
			FINDMINMAX(v0.y,v1.y,v2.y,min,max);
			if(min>boxHalfSize.y || max<-boxHalfSize.y) 
				return false;

			// test in Z-direction 
			FINDMINMAX(v0.z,v1.z,v2.z,min,max);
			if(min>boxHalfSize.z || max<-boxHalfSize.z) 
				return false;

			// Bullet 2: 
			//  test if the box intersects the plane of the triangle 
			//  compute plane equation of triangle: normal*x+d=0 
			normal=cross(e0,e1);

			// -NJMP- (line removed here)
			if(!planeBoxOverlap(normal,v0,boxHalfSize)) 
				return false; // -NJMP-

			return true; // box and triangle overlaps 
		}
		//returns the closest point in the plane
		data3<T> closestPoint(const data3<T> &pos) const
		{
			data3<T> edge0 = vertices[1] - vertices[0];
			data3<T> edge1 = vertices[2] - vertices[0];
			data3<T> v0 = vertices[0] - pos;

			T a = dot(edge0,edge0 );
			T b = dot(edge0,edge1 );
			T c = dot(edge1,edge1 );
			T d = dot(edge0,v0 );
			T e = dot(edge1,v0 );

			T det = a*c - b*b;
			T s = b*e - c*d;
			T t = b*d - a*e;

			if ( s + t < det )
			{
				if ( s < (T)(0.0) )
				{
					if ( t < (T)(0.0) )
					{
						if ( d < (T)(0.0) )
						{
							s = CLAMP( -d/a, (T)(0.0), (T)(1.0) );
							t = (T)(0.0);
						}
						else
						{
							s = (T)(0.0);
							t = CLAMP( -e/c, (T)(0.0), (T)(1.0) );
						}
					}
					else
					{
						s = (T)(0.0);
						t = CLAMP( -e/c, (T)(0.0), (T)(1.0) );
					}
				}
				else if ( t < (T)(0.0) )
				{
					s = CLAMP( -d/a, (T)(0.0), (T)(1.0) );
					t = (T)(0.0);
				}
				else
				{
					T invDet = (T)(1.0) / det;
					s *= invDet;
					t *= invDet;
				}
			}
			else
			{
				if ( s < (T)(0.0) )
				{
					T tmp0 = b+d;
					T tmp1 = c+e;
					if ( tmp1 > tmp0 )
					{
						T numer = tmp1 - tmp0;
						T denom = a-2*b+c;
						s = CLAMP( numer/denom, (T)(0.0), (T)(1.0) );
						t = 1-s;
					}
					else
					{
						t = CLAMP( -e/c, (T)(0.0), (T)(1.0) );
						s = (T)(0.0);
					}
				}
				else if ( t < (T)(0.0) )
				{
					if ( a+d > b+e )
					{
						T numer = c+e-b-d;
						T denom = a-2*b+c;
						s = CLAMP( numer/denom, (T)(0.0), (T)(1.0) );
						t = 1-s;
					}
					else
					{
						s = CLAMP( -e/c, (T)(0.0), (T)(1.0) );
						t = (T)(0.0);
					}
				}
				else
				{
					T numer = c+e-b-d;
					T denom = a-2*b+c;
					s = CLAMP( numer/denom, (T)(0.0), (T)(1.0) );
					t = (T)(1.0) - s;
				}
			}

			return vertices[0] + edge0 * s + edge1 * t;
		}
		//distance from point to line
		T distance(const data3<T> &pos)
		{
			//returns the closest point in the plane
			data3<T> cPoint=closestPoint(pos);

			//return length
			return (T)length(cPoint-pos);
		}
		//computes bounding box
		box<T> getBoundingBox(void) const
		{
			//compute the bounds
			box<T> bounds(
				MIN(vertices[0].x,MIN(vertices[1].x,vertices[2].x)),
				MIN(vertices[0].y,MIN(vertices[1].y,vertices[2].y)),
				MIN(vertices[0].z,MIN(vertices[1].z,vertices[2].z)),
				MAX(vertices[0].x,MAX(vertices[1].x,vertices[2].x)),
				MAX(vertices[0].y,MAX(vertices[1].y,vertices[2].y)),
				MAX(vertices[0].z,MAX(vertices[1].z,vertices[2].z)));
			return bounds;
		}

};
typedef triangle<int> intTriangle;
typedef triangle<float> floatTriangle;
typedef triangle<double> doubleTriangle;

template <typename T>
class face
{
	public:
		triangle<T> vertices;
		triangle<T> textures;
		triangle<T> normals;
		triangle<T> colors;
		unsigned int materialID;	

		face(void)
		{

		}

		face(const data3<T> &vertex1,const data3<T> &vertex2,const data3<T> &vertex3)
		{
			vertices.vertices[0]=vertex1;
			vertices.vertices[1]=vertex2;
			vertices.vertices[2]=vertex3;
		}

		face(const data3<T> &vertex1,const data3<T> &vertex2,const data3<T> &vertex3,
			const data3<T> &texture1,const data3<T> &texture2,const data3<T> &texture3,
			const data3<T> &normal1,const data3<T> &normal2,const data3<T> &normal3,
			const data3<T> &color1,const data3<T> &color2,const data3<T> &color3)
		{
			vertices.vertices[0]=vertex1;
			vertices.vertices[1]=vertex2;
			vertices.vertices[2]=vertex3;

			textures.vertices[0]=texture1;
			textures.vertices[1]=texture2;
			textures.vertices[2]=texture3;

			normals.vertices[0]=normal1;
			normals.vertices[1]=normal2;
			normals.vertices[2]=normal3;

			colors.vertices[0]=color1;
			colors.vertices[1]=color2;
			colors.vertices[2]=color3;
		}

		//clears the data
		void clear(void)
		{
			vertices.clear();
			textures.clear();
			normals.clear();
			colors.clear();
		}

		const bool operator==(const face<T> &data) const
		{
			return (materialID == data.materialID && vertices==data.vertices && textures == data.textures && normals==data.normals && colors==data.colors);
		}
		
		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const face& data)
		{
			out<<"{"<<data.vertices<<","<<data.textures<<","<<data.normals<<","<<data.colors<<","<<data.materialID<<"}";
			return out;
		}

		//checks if the triange intersects with the ray
		bool doesIntersect(const rayData<T> &ray) const
		{
			return vertices.doesIntersect(ray);
		}

		
		//returns barycentric coordinate
		data3<T> getBarycentric(const data3<T> &pos) const
		{
			return vertices.getBarycentric(pos);
		}

		//get returns the intersection distance
		T getIntersectDistance(const rayData<T> &ray) const
		{
			return vertices.getIntersectDistance(ray);
		}

		//checks if the box intersects the face
		bool doesIntersect(const box<T> &boxData) const
		{
			return vertices.doesIntersect(boxData);
		}
		//returns the closest point in the plane
		data3<T> closestPoint(const data3<T> &pos) const
		{
			return vertices.closestPoint(pos);
		}
		//distance from point to line
		T distance(const data3<T> &pos)
		{
			return vertices.distance(pos);
		}
		//computes bounding box
		box<T> getBoundingBox(void) const
		{
			return vertices.getBoundingBox();
		}



};
typedef face<int> intFace;
typedef face<float> floatFace;
typedef face<double> doubleFace;


template <typename T>
class sphere
{
	public:
		//center of sphere
		data3<T> center;

		//radius of the sphere
		T radius;

		uchar3 diffuseColor;
		unsigned int id;

		sphere(void)
		{

		}

		sphere(const data3<T> &center,const T &radius)
		{
			this->center=center;
			this->radius=radius;
		}

		//clears the data
		void clear(void)
		{
			//center of sphere
			center.clear();
			//radius of the sphere
			radius=(T)0;
			diffuseColor.clear();
			id=0;
		}

		//smaller than operator
		bool operator<(const sphere<T> &data) const
		{
			return false;
		}

		// data equality
		const bool operator==(const sphere<T> &data) const
		{
			return (center==data.center && radius == data.radius);
		}

		//distance from point to sphere
		T distance(const data3<T> &pos)
		{
			//return length
			return MIN((T)length(center-pos)-radius,(T)0);
		}

		int planeBoxOverlap(data3<T> &normal, data3<T> &vert, data3<T> &maxbox) // -NJMP-
		{
			return false;
		}


		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const sphere& data)
		{
			out<<"{"<<data.center<<","<<data.radius<<"}";
			return out;
		}

		//checks if the triange intersects with the ray
		bool doesIntersect(const rayData<T> &ray) const
		{
			T a,b,c;
			T temp;
			//a=dir.x*dir.x+dir.y*dir.y+dir.z*dir.z;
			//b=2 * (dir.x * (pos.x - center.x) + dir.y * (pos.y - center.y) + dir.z * (pos.z - center.z));
			//c=(pos.x - center.x)*(pos.x - center.x)+(pos.y - center.y)*(pos.y - center.y)+(pos.z - center.z)*(pos.z - center.z) - radius*radius;

			a=dot(ray.direction,ray.direction);
			b=2*dot(ray.direction,(ray.position-center));
			c=dot((ray.position-center),(ray.position-center))-(radius*radius);

			temp=(b*b)-(4*a*c);
			
			if(temp<0)	//doesn't intersect
				return false;

			if(((-b-(T)sqrt(temp))/(2*a))<=0 && ((-b+(T)sqrt(temp))/(2*a))<=0)
				return false;
			return true;;	//intersects
			
		}

		//get returns the intersection distance
		T getIntersectDistance(const rayData<T> &ray) const
		{
			//float a,b,c;
			//float temp;

			////a=dir.x*dir.x+dir.y*dir.y+dir.z*dir.z;
			////b=2 * (dir.x * (pos.x - center.x) + dir.y * (pos.y - center.y) + dir.z * (pos.z - center.z));
			////c=(pos.x - center.x)*(pos.x - center.x)+(pos.y - center.y)*(pos.y - center.y)+(pos.z - center.z)*(pos.z - center.z) - radius*radius;

			//a=dot(ray.direction,ray.direction);
			//b=2*dot(ray.direction,(ray.position-center));
			//c=dot((ray.position-center),(ray.position-center))-(radius*radius);

			//temp= (b*b)-(4*a*c);

			//if(temp<0)	//doesn't intersect
			//	return (T)-1;

			//float t0,t1;
			//t0=(-b-(float)sqrt(temp))/(2*a);
			//t1=(-b+(float)sqrt(temp))/(2*a);

			//if(t0<=0 && t1<=0)
			//	return (T)-1;;

			//return t0;

			double a,b,c;
			double temp;

			double3 rayDirectionD(ray.direction.x,ray.direction.y,ray.direction.z);
			double3 rayPositionD(ray.position.x,ray.position.y,ray.position.z);
			double3 centerD(center.x,center.y,center.z);
			double radiusD=radius;


			a=dot(rayDirectionD,rayDirectionD);
			b=2*dot(rayDirectionD,(rayPositionD-centerD));
			c=dot((rayPositionD-centerD),(rayPositionD-centerD))-(radiusD*radiusD);

			temp= (b*b)-(4*a*c);

			if(temp<0)	//doesn't intersect
				return (T)-1;

			double t0,t1;
			//t0=(-b-sqrt(temp))/(2*a);
			t1=(-b+sqrt(temp))/(2*a);
			if(t1<=0)
				return (T)-1;

			//since c=t0*t1...look in numerical analysis book for precision related stuff
			t0=c/t1;

			if(t0<=0)
				return (T)-1;

			return (T) t0;
		}

		//computes bounding box
		box<T> getBoundingBox(void) const
		{
			//compute the bounds
			box<T> bounds(center.x-radius,center.y-radius,center.z-radius,center.x+radius,center.y+radius,center.z+radius);
			return bounds;
		}

		//checks if the sphere intersects with the box
		bool doesIntersect(const box<T> &boxData) const
		{

			T s, d =(T) 0; 
			//find the square of the distance
			//from the sphere to the box
			for(unsigned int i=0; i<3 ; i++ ) 
			{ 

				if( center.xyz[i] < boxData.min[i] )
				{ 


					s = center.xyz[i] - boxData.min[i];
					d += s*s; 
				} 

				else if( center.xyz[i] > boxData.max[i] )
				{ 

					s = center.xyz[i] - boxData.max[i];
					d += s*s; 
				} 
			}
			return d <= radius*radius; 
		}

};
typedef sphere<int> intSphere;
typedef sphere<float> floatSphere;
typedef sphere<double> doubleSphere;

template <typename T>
class frustum
{
		//extract a plane from a given matrix and row id
		void extractPlane(data4<T> &plane, data4x4<T> &mat, int row)
		{
			int scale = (row < 0) ? -1 : 1;
			row = abs(row) - 1;

			// calculate plane coefficients from the matrix
			plane.x = mat[3] + scale * mat[row];
			plane.y = mat[7] + scale * mat[row + 4];
			plane.z = mat[11] + scale * mat[row + 8];
			plane.w = mat[15] + scale * mat[row + 12];

			// normalize the plane by using only first 3 coefficients
			T length = (T) sqrt(plane.x * plane.x + plane.y * plane.y + plane.z * plane.z);
			if(length)
				plane/=length;

		}

	public:

		//clipping planes
		data4<T> planes[6];

		//clear the data
		void clear(void)
		{
			planes[0].clear();
			planes[1].clear();
			planes[2].clear();
			planes[3].clear();
			planes[4].clear();
			planes[5].clear();
		}

		//calculate frustum
		void calculateFrustum(const data4x4<T> &modelViewProj)
		{
			//calculate modelview projection matrix
			data4x4<T> modelViewProjTranspose=transpose(modelViewProj);
				
			// extract each plane
			//left
			extractPlane(planes[0], modelViewProjTranspose, 1);
			
			//right
			extractPlane(planes[1], modelViewProjTranspose, -1);
			
			//bottom
			extractPlane(planes[2], modelViewProjTranspose, 2);
			
			//top
			extractPlane(planes[3], modelViewProjTranspose, -2);
			
			//near
			extractPlane(planes[4], modelViewProjTranspose, 3);
			
			//far
			extractPlane(planes[5], modelViewProjTranspose, -3);
		}

		//returns true if point is inside the camera frustum
		bool doesIntersect(const data3<T> &vec) const
		{
			int i;
			for ( i = 0; i < 6; ++i)
			{
				if(( planes[i].x * vec.x + planes[i].y *vec.y +	planes[i].z * vec.z + planes[i].w) < 0 )
					return false;
			}
			return true;
		}

		//returns true if the sphere is inside the frustum
		bool doesIntersect(const sphere<T> &sphere) const
		{
			if (planes[0].x * sphere.center.x +	planes[0].y * sphere.center.y +	planes[0].z * sphere.center.z +	planes[0].w	 < -sphere.radius)
				return false;
			if (planes[1].x * sphere.center.x +	planes[1].y * sphere.center.y +	planes[1].z * sphere.center.z +	planes[1].w	 < -sphere.radius)
				return false;
			if (planes[2].x * sphere.center.x +	planes[2].y * sphere.center.y +	planes[2].z * sphere.center.z +	planes[2].w	 < -sphere.radius)
				return false;
			if (planes[3].x * sphere.center.x +	planes[3].y * sphere.center.y +	planes[3].z * sphere.center.z +	planes[3].w	 < -sphere.radius)
				return false;
			if (planes[4].x * sphere.center.x +	planes[4].y * sphere.center.y +	planes[4].z * sphere.center.z +	planes[4].w	 < -sphere.radius)
				return false;
			if (planes[5].x * sphere.center.x +	planes[5].y * sphere.center.y +	planes[5].z * sphere.center.z +	planes[5].w	 < -sphere.radius)
				return false;
			return true;
		}

		//returns true if the sphere is inside the frustum
		bool doesIntersect(const data3<T> &center, const T &radius) const
		{
			if (planes[0].x * center.x +	planes[0].y * center.y +	planes[0].z * center.z +	planes[0].w	 < -radius)
				return false;
			if (planes[1].x * center.x +	planes[1].y * center.y +	planes[1].z * center.z +	planes[1].w	 < -radius)
				return false;
			if (planes[2].x * center.x +	planes[2].y * center.y +	planes[2].z * center.z +	planes[2].w	 < -radius)
				return false;
			if (planes[3].x * center.x +	planes[3].y * center.y +	planes[3].z * center.z +	planes[3].w	 < -radius)
				return false;
			if (planes[4].x * center.x +	planes[4].y * center.y +	planes[4].z * center.z +	planes[4].w	 < -radius)
				return false;
			if (planes[5].x * center.x +	planes[5].y * center.y +	planes[5].z * center.z +	planes[5].w	 < -radius)
				return false;
			return true;
		}

		//returns true if one of the edge is inside the frustum
		//becareful while using this, it has some limitation
		bool doesIntersect(const data3<T> *vec) const
		{
			int i,j;
			bool point;
			for(j=0;j<8;j++)
			{
				point=true;
				//if for all the plane the number in >0 then the point is inside
				for ( i = 0; i < 6; i++)
				{
					if(( planes[i].x * vec[j].x + planes[i].y *
						vec[j].y +	planes[i].z * vec[j].z + planes[i].w) < 0 )
					{
						point=false;
						break;
					}
				}
				if(point)
					return true;
			}
			return false;
		}

		//returns true if the cube is inside the frustum
		bool doesIntersectCube(const sphere<T> &sphere) const
		{
			data3<T> *ptr=&sphere.center;
			T &distance=sphere.radius;
			for(int i = 0; i < 6; i++ )
			{
				if(planes[i].x * (ptr->x - distance) + planes[i].y * (ptr->y - distance) + planes[i].z * (ptr->z - distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x + distance) + planes[i].y * (ptr->y - distance) + planes[i].z * (ptr->z - distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x - distance) + planes[i].y * (ptr->y + distance) + planes[i].z * (ptr->z - distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x + distance) + planes[i].y * (ptr->y + distance) + planes[i].z * (ptr->z - distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x - distance) + planes[i].y * (ptr->y - distance) + planes[i].z * (ptr->z + distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x + distance) + planes[i].y * (ptr->y - distance) + planes[i].z * (ptr->z + distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x - distance) + planes[i].y * (ptr->y + distance) + planes[i].z * (ptr->z + distance) + planes[i].w > 0)
					continue;
				if(planes[i].x * (ptr->x + distance) + planes[i].y * (ptr->y + distance) + planes[i].z * (ptr->z + distance) + planes[i].w > 0)
					continue;
				// if we get here, it isn't in the frustum
				return false;
			}			
			return true;
		}
};
typedef frustum<int> intFrustum;
typedef frustum<float> floatFrustum;
typedef frustum<double> doubleFrustum;

template <typename T>
class quadrilateral
{
	public:
		T minX,minY,maxX,maxY;

		quadrilateral(void)	{}

		quadrilateral(const T &minX,const T &minY,const T &maxX,const T &maxY)
		{
			this->minX=minX;
			this->minY=minY;
			this->maxX=maxX;
			this->maxY=maxY;
		}

		//clears the data
		void clear(void)
		{
			minX=minY=maxX=maxY=0;
		}

		bool doesIntersect(const float2 &point)
		{
			return minX<=point.x && point.x<maxX && minY<=point.y && point.y<maxY;
		}

		//closest distance between quad and the position
		T distance(const data2<T> &pos)
		{
			//create 4 lines
			float2Line a(float2S(minX,minY),float2S(maxX,minY));
			float2Line b(float2S(maxX,minY),float2S(maxX,maxY));
			float2Line c(float2S(maxX,maxY),float2S(minX,maxY));
			float2Line d(float2S(minX,maxY),float2S(minX,minY));

			T dist=a.distance(pos);
			dist=MIN(dist,b.distance(pos));
			dist=MIN(dist,c.distance(pos));
			dist=MIN(dist,d.distance(pos));

			return dist;
		}


		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const quadrilateral& data)
		{
			out<<"<"<<data.minX<<", "<<data.minY<<", "<<data.maxX<<", "<<data.maxY<<">";
			return out;
		}
};
typedef quadrilateral<int> intQuad;
typedef quadrilateral<float> floatQuad;
typedef quadrilateral<double> doubleQuad;

template <typename T>
class box
{
	public:
		

		union
		{
			struct
			{
				T minX,minY,minZ,maxX,maxY,maxZ;
			};
					
			struct
			{
				T min[3];
				T max[3];
			};
		};

		box(void){}

		box(const T &minX,const T &minY,const T &minZ,const T &maxX,const T &maxY,const T &maxZ)
		{
			this->minX=minX;
			this->minY=minY;
			this->minZ=minZ;
			this->maxX=maxX;
			this->maxY=maxY;
			this->maxZ=maxZ;
		}

		//clears the data
		void clear(void)
		{
			minX=minY=minZ=maxX=maxY=maxZ=0;
		}

		//checks if the point is insie the box
		bool doesIntersect(const float3 &point)
		{
			return minX<=point.x && point.x<maxX && minY<=point.y && point.y<maxY && minZ<=point.z && point.z<maxZ;
		}

		//checks if the point is insie the box
		bool doesIntersect(const floatPhotonData &pData)
		{
			return minX<=pData.position.x && pData.position.x<maxX && minY<=pData.position.y && pData.position.y<maxY && minZ<=pData.position.z && pData.position.z<maxZ;
		}

		//checks if the sphere intersects with the box
		bool doesIntersect(const sphere<T> &sph) const
		{
			return sph.doesIntersect(*this);
		}

		//retruns the half Size of the box
		data3<T> getHalfSize(void) const
		{
			return data3<T>((T)((maxX-minX)*.5),(T)((maxY-minY)*.5),(T)((maxZ-minZ)*.5));
		}

		//retruns the center of the box
		data3<T> getCenter(void) const
		{
			return data3<T>((T)((maxX+minX)*.5),(T)((maxY+minY)*.5),(T)((maxZ+minZ)*.5));
		}

		//distance from point to plane
		T distance(const data3<T> &pos) const
		{
			//make 6 planes
			plane<T> a(data3<T>(minX,minY,minZ),data3<T>(minX,maxY,minZ),data3<T>(minX,minY,maxZ));
			plane<T> b(data3<T>(minX,minY,maxZ),data3<T>(minX,maxY,maxZ),data3<T>(maxX,maxY,maxZ));
			plane<T> c(data3<T>(maxX,minY,minZ),data3<T>(maxX,minY,maxZ),data3<T>(maxX,maxY,minZ));
			plane<T> d(data3<T>(minX,minY,minZ),data3<T>(maxX,minY,minZ),data3<T>(maxX,maxY,minZ));
			plane<T> e(data3<T>(minX,minY,minZ),data3<T>(minX,minY,maxZ),data3<T>(maxX,minY,maxZ));
			plane<T> f(data3<T>(minX,maxY,minZ),data3<T>(maxX,maxY,maxZ),data3<T>(minX,maxY,maxZ));

			//get smallest distance
			T dist=a.distance(pos);
			dist=MIN(dist,b.distance(pos));
			dist=MIN(dist,c.distance(pos));
			dist=MIN(dist,d.distance(pos));
			dist=MIN(dist,e.distance(pos));
			dist=MIN(dist,f.distance(pos));

			return dist;
		}

		//======================== X-tests ========================
		#define AXISTEST_X01(a, b, fa, fb)           \
		   p0 = a*v0.y - b*v0.z;                   \
		   p2 = a*v2.y - b*v2.z;                   \
		   if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
		   rad = fa * boxHalfSize.y + fb * boxHalfSize.z;   \
		   if (min>rad || max<-rad) return false;

		#define AXISTEST_X2(a, b, fa, fb)            \
		   p0 = a*v0.y - b*v0.z;                  \
		   p1 = a*v1.y - b*v1.z;                   \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.y + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		//======================== Y-tests ========================
		#define AXISTEST_Y02(a, b, fa, fb)           \
		   p0 = -a*v0.x + b*v0.z;                  \
		   p2 = -a*v2.x + b*v2.z;                        \
		   if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		#define AXISTEST_Y1(a, b, fa, fb)            \
		   p0 = -a*v0.x + b*v0.z;                  \
		   p1 = -a*v1.x + b*v1.z;                     \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.z;   \
		   if(min>rad || max<-rad) return false;

		//======================== Z-tests ========================
		#define AXISTEST_Z12(a, b, fa, fb)           \
		   p1 = a*v1.x - b*v1.y;                  \
		   p2 = a*v2.x - b*v2.y;                   \
		   if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.y;   \
		   if(min>rad || max<-rad) return false;


		#define AXISTEST_Z0(a, b, fa, fb)            \
		   p0 = a*v0.x - b*v0.y;             \
		   p1 = a*v1.x - b*v1.y;                  \
		   if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
		   rad = fa * boxHalfSize.x + fb * boxHalfSize.y;   \
		   if(min>rad || max<-rad) return false;

		#define FINDMINMAX(x0,x1,x2,min,max) \
		   min = max = x0;   \
		   if(x1<min) min=x1;\
		   if(x1>max) max=x1;\
		   if(x2<min) min=x2;\
		   if(x2>max) max=x2;

		//checks if the box intersects with the triangle
		bool doesIntersect(const triangle<T> &triangleData)
		{
			data3<T> boxHalfSize=getHalfSize();
			data3<T> boxCenter=getCenter();

			//use separating axis theorem to test overlap between triangle and box 
			//    need to test for overlap in these directions: 
			//    1) the {x,y,z}-directions (actually, since we use the AABB of the triangle 
			//       we do not even need to test these) 
			//    2) normal of the triangle 
			//    3) crossproduct(edge from tri, {x,y,z}-directin) 
			//       this gives 3x3=9 more tests 
			data3<T> v0,v1,v2;
			T min,max,p0,p1,p2,rad,fex,fey,fez;      // -NJMP- "d" local variable removed
			data3<T> normal,e0,e1,e2;

			// This is the fastest branch on Sun 
			// move everything so that the boxCenter is in (0,0,0) 
			v0=triangleData.vertices[0]-boxCenter;
			v1=triangleData.vertices[1]-boxCenter;
			v2=triangleData.vertices[2]-boxCenter;

			// compute triangle edges 
			e0=v1-v0;      // tri edge 0 
			e1=v2-v1;      // tri edge 1 
			e2=v0-v2;      // tri edge 2 

			// Bullet 3:  
			//  test the 9 tests first (this was faster) 
			fex = abs(e0.x);
			fey = abs(e0.y);
			fez = abs(e0.z);
			AXISTEST_X01(e0.z, e0.y, fez, fey);
			AXISTEST_Y02(e0.z, e0.x, fez, fex);
			AXISTEST_Z12(e0.y, e0.x, fey, fex);

			fex = abs(e1.x);
			fey = abs(e1.y);
			fez = abs(e1.z);
			AXISTEST_X01(e1.z, e1.y, fez, fey);
			AXISTEST_Y02(e1.z, e1.x, fez, fex);
			AXISTEST_Z0(e1.y, e1.x, fey, fex);

			fex = abs(e2.x);
			fey = abs(e2.y);
			fez = abs(e2.z);
			AXISTEST_X2(e2.z, e2.y, fez, fey);
			AXISTEST_Y1(e2.z, e2.x, fez, fex);
			AXISTEST_Z12(e2.y, e2.x, fey, fex);

			// Bullet 1: 
			//  first test overlap in the {x,y,z}-directions 
			//  find min, max of the triangle each direction, and test for overlap in 
			//  that direction -- this is equivalent to testing a minimal AABB around 
			//  the triangle against the AABB 

			// test in X-direction 
			FINDMINMAX(v0.x,v1.x,v2.x,min,max);
			if(min>boxHalfSize.x || max<-boxHalfSize.x)
				return false;

			// test in Y-direction 
			FINDMINMAX(v0.y,v1.y,v2.y,min,max);
			if(min>boxHalfSize.y || max<-boxHalfSize.y) 
				return false;

			// test in Z-direction 
			FINDMINMAX(v0.z,v1.z,v2.z,min,max);
			if(min>boxHalfSize.z || max<-boxHalfSize.z) 
				return false;

			// Bullet 2: 
			//  test if the box intersects the plane of the triangle 
			//  compute plane equation of triangle: normal*x+d=0 
			normal=cross(e0,e1);

			// -NJMP- (line removed here)
			if(!planeBoxOverlap(normal,v0,boxHalfSize)) 
				return false; // -NJMP-

			return true; // box and triangle overlaps 
		}

		//checks if the box intersects with face
		bool doesIntersect(const face<T> &fData)
		{
			return doesIntersect(fData.vertices);
		}
		//checks if the box intersects with the triangle
		//Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley
		bool doesIntersect(const rayData<T> &ray) const
		{
			data3<T> bbox[2];
			
			bbox[0].x=minX;
			bbox[0].y=minY;
			bbox[0].z=minZ;

			bbox[1].x=maxX;
			bbox[1].y=maxY;
			bbox[1].z=maxZ;
			
			data3<T> invDir((T)(1.0/ray.direction.x),(T)(1.0/ray.direction.y),(T)(1.0/ray.direction.z));
			data3<unsigned int> sign((unsigned int)(invDir.x<0),(unsigned int)(invDir.y<0),(unsigned int)(invDir.z<0));

			T tmin = (bbox[sign.x].x - ray.position.x) * invDir.x;
			T tmax = (bbox[1-sign.x].x - ray.position.x) * invDir.x;

			T tymin = (bbox[sign.y].y - ray.position.y) * invDir.y;
			T tymax = (bbox[1-sign.y].y - ray.position.y) * invDir.y;

			if ( (tmin > tymax) || (tymin > tmax) )
				return false;
			if (tymin > tmin)
				tmin = tymin;
			if (tymax < tmax)
				tmax = tymax;

			T tzmin = (bbox[sign.z].z - ray.position.z) * invDir.z;
			T tzmax = (bbox[1-sign.z].z - ray.position.z) * invDir.z;

			if ( (tmin > tzmax) || (tzmin > tmax) )
				return false;
			if (tzmin > tmin)
				tmin = tzmin;
			if (tzmax < tmax)
				tmax = tzmax;
			if(tmin<0 && tmax<0)
				return false;
			return true;
		}

		//get returns the intersection distance
		T getIntersectDistance(const rayData<T> &ray) const
		{
			data3<T> bbox[2];
			
			bbox[0].x=minX;
			bbox[0].y=minY;
			bbox[0].z=minZ;

			bbox[1].x=maxX;
			bbox[1].y=maxY;
			bbox[1].z=maxZ;
			
			data3<T> invDir((T)(1.0/ray.direction.x),(T)(1.0/ray.direction.y),(T)(1.0/ray.direction.z));
			data3<unsigned int> sign((unsigned int)(invDir.x<0),(unsigned int)(invDir.y<0),(unsigned int)(invDir.z<0));

			T tmin = (bbox[sign.x].x - ray.position.x) * invDir.x;
			T tmax = (bbox[1-sign.x].x - ray.position.x) * invDir.x;

			T tymin = (bbox[sign.y].y - ray.position.y) * invDir.y;
			T tymax = (bbox[1-sign.y].y - ray.position.y) * invDir.y;

			if ( (tmin > tymax) || (tymin > tmax) )
				return -1;
			if (tymin > tmin)
				tmin = tymin;
			if (tymax < tmax)
				tmax = tymax;

			T tzmin = (bbox[sign.z].z - ray.position.z) * invDir.z;
			T tzmax = (bbox[1-sign.z].z - ray.position.z) * invDir.z;

			if ( (tmin > tzmax) || (tzmin > tmax) )
				return -1;
			if (tzmin > tmin)
				tmin = tzmin;
			if (tzmax < tmax)
				tmax = tzmax;

			if(tmin<0)
			{
				if(tmax<0)
					return -1;
				return tmax;
			}
			return tmin;
		}

		int planeBoxOverlap(data3<T> &normal, data3<T> &vert, data3<T> &maxbox) // -NJMP-
		{
			data3<T> vmin,vmax;

			if(normal.x>(T)0.0)
			{
				vmin.x=-maxbox.x - vert.x; // -NJMP-
				vmax.x= maxbox.x - vert.x; // -NJMP-
			}
			else
			{
				vmin.x= maxbox.x - vert.x; // -NJMP-
				vmax.x=-maxbox.x - vert.x; // -NJMP-
			}

			if(normal.y>(T)0.0)
			{
				vmin.y=-maxbox.y - vert.y; // -NJMP-
				vmax.y= maxbox.y - vert.y; // -NJMP-
			}
			else
			{
				vmin.y= maxbox.y - vert.y; // -NJMP-
				vmax.y=-maxbox.y - vert.y; // -NJMP-
			}

			if(normal.z>(T)0.0)
			{
				vmin.z=-maxbox.z - vert.z; // -NJMP-
				vmax.z= maxbox.z - vert.z; // -NJMP-
			}
			else
			{
				vmin.z= maxbox.z - vert.z; // -NJMP-
				vmax.z=-maxbox.z - vert.z; // -NJMP-
			}

			if(dot(normal,vmin)>(T)0.0) 
				return false;  // -NJMP-
			if(dot(normal,vmax)>=(T)0.0) 
				return true; // -NJMP-

			return false;
		}

		//prints to the output stream
		friend std::ostream& operator<<(std::ostream &out, const box& data)
		{
			out<<"<"<<data.minX<<", "<<data.minY<<", "<<data.minZ<<", "<<data.maxX<<", "<<data.maxY<<", "<<data.maxZ<<">";
			return out;
		}
};
typedef box<int> intBox;
typedef box<float> floatBox;
typedef box<double> doubleBox;


//data structure to hold setup data
class setupData
{
	public:
		char *name;
		char *value;

		setupData(void):name(NULL),value(NULL)
		{
		}

		setupData(char *&namePtr,char *&valuePtr)
		{
			name=NULL;
			value=NULL;
			
			if(namePtr && strlen(namePtr))
			{
				name=new char[strlen(namePtr)+1];
				strcpy(name,namePtr);
			}

			if(valuePtr && strlen(valuePtr))
			{
				value=new char[strlen(valuePtr)+1];
				strcpy(value,valuePtr);
			}
		}

		setupData(const setupData &sData)
		{
			name=NULL;
			value=NULL;
			
			if(sData.name && strlen(sData.name))
			{
				name=new char[strlen(sData.name)+1];
				strcpy(name,sData.name);
			}

			if(sData.value && strlen(sData.value))
			{
				value=new char[strlen(sData.value)+1];
				strcpy(value,sData.value);
			}
		}

		//assignment operator
		const setupData& operator=(const setupData &sData)
		{
			if(this!=((const setupData *)&sData))
			{
				clear();
			
				if(sData.name && strlen(sData.name))
				{
					name=new char[strlen(sData.name)+1];
					strcpy(name,sData.name);
				}

				if(sData.value && strlen(sData.value))
				{
					value=new char[strlen(sData.value)+1];
					strcpy(value,sData.value);
				}
			}
			return *this;
		}

		void clear(void)
		{
			delete[] name;
			delete[] value;
			name=NULL;
			value=NULL;
		}

		~setupData()
		{
			clear();
		}
};


#endif
